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NOMOGRAMS FOR THE SOLUTION OF THE 
SOUND-RANGING PROBLEM IN A PLANE 
By G. C. CURTIS (King’s College, Newcastle-upon-T yne) 
[Received 3 April 1952; in revised form 6 February 1953] 


SUMMARY 

Nomograms are devised for locating a source of sound in the plane defined by three 
collinear receivers and the source, given the times of arrival of the sound at the 
receivers. This is formally equivalent to the construction of a circle to touch three 
given circles having collinear centres. 

Though the time of emission of the sound is not usually of interest, its determina- 
tion forms a useful intermediate stage in the calculations. Thus three unknowns 
must be determined, but it is shown that only two nomograms are required: the first 
giving the time of emission and one coordinate of the source by a single setting of 
the index-line, and the second giving the remaining coordinate. 

The nomograms described possess special advantages when the solution is over- 
determined by the use of data from more than three receivers, since each additional 
receiver may be dealt with by adding an additional scale to the first nomogram. 
The mean solution is then obtained by setting the index along the best straight line 
through the points denoted by the entry values of all the observed distance differ- 


ences. The second nomogram is unaffected. 


1. Introduction 

Tue general problem of sound-ranging is the calculation of the position of 
the source of a sudden sound from observations of the difference in time of 
arrival of the sound at receivers in known positions, the actual time of 
emission being unknown. The velocity of sound is required, to convert the 
observed time-differences into distance-differences: the subject of this 
article is the next step in the calculation, in which the coordinates of the 
source are found. 

In the case to be considered, the source and the receiver lie in a fixed plane. 
If in a scale plot (Fig. 1) we draw circles about the receivers having radii 
exhibiting the observed distance-differences, and construct an additional 
circle to touch them, the source of sound will be found at its centre. Such 


a graphical solution would be of limited accuracy, more especially if data 


from more than the minimum of three receivers were to be effectively 
utilized. 

Alternatively, with each pair of receivers as foci a family of hyperbolae 
could be drawn, each labelled with a value of the distance difference. The 
source would be found at the intersection of the appopriate hyperbolae, one 
from each family. The labour of constructing such charts would increase 
rapidly with the accuracy required, but might be justified if usage were 
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130 G. C. CURTIS 

heavy (e.g. in the analogous radar case of the Decca Navigator). These 

two methods have the merit, however, of being usable for non-collinear 

receivers, and either might furnish a first approximation for iterative com- 

putation in such a case, for which direct calculation is very unwieldy. 
Collinear receiver-arrays give rise to simpler computation, but it is still 

desirable to find means of lessening the labour involved. 






P(6ource of sound) 





Fic. 1. Location of a source of sound by means of the collinear receivers M, to M,. 


It is understood that charts for the computation are published in 
restricted Service literature. Those to whom these charts are not available 
may find the method to be described useful: it may also prove the more 
practical in certain cases, for example, when the receivers are not equally 
spaced. An electrical analogue computer would be an alternative (1). 


2. Equations to be solved 

In Fig. 1, P is the sound source to be located and M, to M, are five collinear 
receivers (not necessarily equally spaced). As three receivers suffice to locate 
a source, only M/,, M,, and M, will at first be considered. 

Take M, as origin and M, M, M, as the x-axis and let the coordinates of 
M,, M,, and My be (x,,0), (0,0), and (a3, 0) respectively. 

If PM, = r then PM, = r+s, and PM, = r+sg, where s, and s, are the 
observed distance differences. 

The equations to be solved for x and y (the coordinates of P) are then 


erty? = 72, (1) 
(~x—2,)?+y? . (r+s,)?, (2) 
(x—ay)?-+y? = (r+). (3) 


3. Nomogram for x and r 
Subtracting (1) from (2) and (3), we have 
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and solving for —x and —2r we find 
X —2r —]) 
2 od 2 — = ain 1” 
8; %%— 8] Yy— Sy hy) yy | 
. 2 2 2 S * Qe.) i dy 
§3  t3—S3 r3—83 =X3| {83 +X] 
or, considering x and r separately, 
v S 72 — 3? | 
S yA d ang | > 
ft “ee SS ee (6x) 
> ») > 4 »2 
og akg [83 %g—83 
is | y2_g2 Dy 
oi* ~"*1 451 = . | - 0. (6r) 
. dap | ~S__ a8 > 
83 <3 | |%3— 83 stg | 


Expanded, (6x) and (6r) would give formulae for computation of x and r. 
To derive nomograms, write them respectively 


0 x ] | 2r 0 
2 »2 Dy 0 ‘ > “2 ee Dy —_ 

S; Xj Sj 2X and 8} xy 8} 22 0. 
2 »2 oy > 2 »2 ay 

S3 %3—S83 2X3 83 %—83 <X3 


Adding in each ease the first column to the last, and then dividing each row 


throughout by its last member, we have the standard nomographic forms 


0 U l l 2r ] 
) > 8 y2 §? 

s Us Sy © J _ 
> 5 l 0 (7x) and : -_* 61) =0. (Tr) 
“ny Sy 2X} $; et | at, 8; | 

, 2 2 “2 22 

83 4S 83 - - 2 
) yy dy > ) adit 
2X5 Se 2X2 Se od 3 183 ya 3 83 


The first nomogram can now be constructed (Fig. 2). With axes OA and 


OB, plot scales of x, r, s,, and s, such that the graduation x is at the point 


» 
P.(0,x), the graduation r at P.(1,2r), the graduation s, at 


, -_— 
P,| | 2 =| 
yee © Mie re 
24,18, 2%,+8, 
and the graduation s, at 
Ss. x2— 8? 
P 3 3 3 
3197 1 6 r,+8 ' 
ot3 3 atg 83 


Then the condition for P, P, P, to be a straight line is (7x), which is the 
solution of (1), (2), and (3) for 2 in terms of s, and s,. Hence a straight-edge 
laid across the scale-readings of s, and s, will intersect the scale of x at the 
required value. By (7r) it will also intersect the scale of r at the required 
value, Also, (7a) and (7r) will remain valid if we multiply the first column 
of each by K and the second by L, where K and L are any convenient 


constants. 


By suitable choice of K and L the nomogram may be given any 
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desired dimensions. In Fig. 2, OC and OD are measures of K and L 
respectively. 
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Fic. 2, Nomogram for determining 2 and 7 from s, and 83. 


4. Nomogram for y 

Once x and r are known, y can be determined from (1), (2), or (3). The 
nomogram (three parallel straight scales, Fig. 3) is almost self-evident and 
an appropriate determinantal form of (1) will be quoted without detailed 
derivation : 


0 2 1) 
1 }? 1/=0, (8) 
|2 y? 1] 


The first two columns should as before be multiplied by suitable scale- 
factors. 

This nomogram may be replaced by a form of slide-rule in which distance 
along the rule is proportional to the square of the scale-reading. In the 
author’s experience, the slide-rule is more troublesome to construct than 
the nomogram but is much more convenient to use, especially when high 
accuracy is required. In his opinion, this should also be true of all equations 
which can be converted into the form f(x)+-g(y)+-h(z) = constant and also, 
with greater emphasis, of corresponding equations with four variables. (In 
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the last case, the nomogram would normally be of the double alignment or 
set-square index type (2).) 











t scale F Scale 
of ac fy 
Scale 
+10- “lr 1-0 
1:24 
+08 - -0°8 
Se ee 2 os 10-4 4q 
+06- 7 ie 
+0-4- 06 4 0-4 
4 0-44 - 
+024 J 0-2 
‘i - - - - - - - = - - - - ° ° - ° - =* -_-* 


Fic. 3. Nomogram for determining y from « and r. 


5. Incorporation of results from extra receivers 

Suppose results from two extra receivers .M, and M, (Fig. 1) are obtained. 
Combine them with the result from JJ/,, temporarily ignoring M, and 3. 
The retention of MV, is essential if we are to use the old scale of r as explained 
below. ‘Feed’ them into the nomogram plotted from versions of (77) and 
ir) in which s, is replaced by 4, 83 by s;, x, by 24, and 2, by 2, (Fig. 4). 
The 2 and r seales will be unchanged, hence the straight-edge P, P, P,P. 
should also pass through the points P, and P; which define s, and s, respec- 
tively. In practice, then, one would take the best straight line through the 
experimental positions of P,, P;, P,, and P;, and read off x and r from it. 
This pair of values can then be used in determining y from Fig. 3, which 
needs no alteration. 

Fig. 4 can be extended to deal with any number of collinear receivers: 
however, it seems best to choose as the ‘singular’ receiver (such as M, in 
this instance) one which lies near the centre of the array. 

The positions of the straight-edges shown in Figs. 2, 3, and 4 correspond 
to the point P in Fig. 1, the relevant data being 
l, v7, = 2, x, = 3 units 
x = 0-7, y = 0-6 units 
Ls. = 0-67, r+s, = 1-43, r+s, = 2-38 units 


S; 0-88, 8, 0-25, s, = 0-51, 8; = 1-46 units. 
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Fic. 4. Nomogram for determining x and r from 8, 83, 84, 

and s;. The graduations at E and F are, respectively, +1 

and —1 for s,, —1 and +1 for s,, —2 and +2 for s,, and 
-3 and +3 for s;. 


6. Conclusion 

The following advantages are claimed for the nomograms described 
above: 

1. The fact that the complete calculation is done by two nomograms, in 
each of which only one setting of the index-line is required, should lead to 
speed and convenience in the field. (If a constant value may be assumed 
for the velocity of sound, time-differences may be directly used by a suitable 
change of graduation.) 

2. The incorporation of results from more than three receivers is at once 
easy and satisfactory, for the standard method of combining results subject 
to errors is used, that is, the choice of the best straight line through a number 
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of points. It is not claimed that this represents a least-squares solution, but 
computation (unless prohibitively elaborate) is much less satisfactory in 
this respect. 

The defect that the observations from the common receiver M, are unduly 
weighted is a frequent one in sound-ranging, and might here be met by 
constructing additional nomograms based on a different common receiver. 
These extra nomograms would be essential if the observations from M, were 
to be absent or unreliable in particular instances. 
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TECHNIQUES FOR THE SUMMATION OF PRODUCTS 
ON HOLLERITH AND NATIONAL ACCOUNTING 
MACHINES 


By ANDREW YOUNG 
(Mathematical Institute, University of Liverpool) 
[Received 17 July 1952] 


SUMMARY 
Methods of evaluating sums of products applicable to numerical quadrature, 
Fourier analysis and synthesis, the evaluation of serial correlations, and other 
computations with Hollerith punched card equipment, are described briefly. One 
method is substantially new. Also given are the modifications required for making 


the methods applicable for use with National accounting machines, the use of 


which, in this type of work, is entirely new. Emphasis is placed throughout on the 
most economical use of machines and on the incorporation of checks. 

1. Introduction 

THE necessity of evaluating extensive sums of products arises in numerical 
quadrature, solutions by least squares, Fourier analysis and synthesis, the 
evaluation of serial correlations, and other mathematical computations. 
Frequently > a?, > ab, > ac,... are all required, the data a appearing in 
each summation. Several specially constructed single-purpose machines 
have been devised for this work—see e.g. (1)—but it is more desirable that 
general-purpose machines should be used. 

Many mathematicians have access to Hollerith punched card equipment 
or to National accounting machines, and methods for summing products 
on the former are well known, but their adaptation for use with National 
accounting machines does not seem to have been published before. Owing 
to the fact that information is conveyed to the Hollerith machine mechani- 
cally, and to the National machine by hand-setting on a keyboard, there 
are changes in emphasis which affect the relative efficiencies of the various 
methods, but, in principle, the machines are so basically similar that any 
method that can be applied on a Hollerith tabulator can be used on a 
National, subject to the overriding considerations of the machines’ counter- 
apacity and distribution facilities. 

In this paper the well-known punched card methods of ‘progressive digit- 
ing’ and its variants and of ‘binary distribution’ are briefly described. The 
method of ‘binary self-rolling’ is described in more detail since it is sub- 
stantially new, although a method depending on the same principle has 
been given by Greenhalgh (2) for use in Fourier synthesis. Modifications 
of all the methods are given for use with National accounting machines. In 
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this case, the methods admit of combination, and as an example of this, 
a technique for summing triple products is described (see section 6 (5)). 
Allthe methods have been used extensively by the author, and in particular, 
the periodogram analyses of two similar time-series have been performed 
on both types of equipment. The times taken are given in an appendix. 
Throughout, emphasis is laid on the incorporation of adequate checks, 


snd also on economy in the use of cards and machines. 


2. Bases of mass-production summations of products 


The bases of the methods available are exemplified by the sum 
23w-+ 14a+ 10y+5z, (1) 


which may be expressed in the alternative forms 


l0fw+(w+a+y)}4+{z+(z+2)+(z+2+w)+4+ (2+2+w)+(z2+2+w)}, (2) 
l6w+ 8(a+y)+4(wta+z)+2(wt+aet+y)+(w+z), (3) 
and 2{2[2(2w+a+y)+(w+a+z)|+(w+a+y)}+w+z. (4) 


On a desk machine, (1) is suitable and quickest; (2), (3), and (4) lead to the 
methods of ‘progressive digiting’, ‘binary distribution’, and ‘binary self- 
rolling’ respectively. 

A further method, ‘constant multiplier distribution’, which is useful in 
1 limited number of applications, is exemplified by numerical quadrature 
formulae of the type 


6 
. 


w 


f(x) dx a0 tl (fo +- fg) +216(f,+-f5)+27(fo+fy)+272f3}, (5) 


0 
the constant multipliers in this case being 


(41, 216, 27, 272)w/140. 


3. Hollerith processes 

The Hollerith equipment required—a tabulator and a mechanical sorter— 
is described by Hartley (3). 

(a) Progressive digiting 

Hartley describes this method and illustrates it in connexion with the 
calculation of serial correlations. In his card layout for the latter the 
iaximum amount of data is not put on to the cards. This does not matter 
when only low-order correlations are required, but if all correlations up to 
i high order are required, as may occur in the preparation of an extensive 
correllogram, the alternative layout giving several times the number of 











138 A. YOUNG 


correlations per card, which is given in Appendix I, is useful. The main 
variants of the method in use are: 

(i) Formation of sums of products six at a time (or twelve at a time if 
the sums are small enough to allow two to be formed in each counter), 
and summary punching the successive digital subtotals on to new 
cards. In this method, by suitable plugging, the summary cards can 
be prepared so that all the cards can be tabulated in one run to give 
the final summations. 

(ii) Formation of sums three at a time (or six if two are formed in each 
counter), and rolling progressive digital subtotals into the remaining 
counters. The unit, tens,... totals are then combined on desk 
machines. 

(iii 


— 


As (i) but without summary punching, the totals being obtaiy — on 
desk machines. 

For the occasional user who has not direct access to a Hollerith installation, 
(iii) is usually most economical. In other cases it is easy to estimate which 
of (i) or (ii) is the better method in a particular computation by considering 
the number of control changes, card passages, desk operations required, 
and also the card consumption. 


(b) Binary distributiont 

In this method, cards are tabulated in such a way that the multiplicands 
are distributed into any of one of four counters according to the value of 
the multiplier as follows: 








Value of multiplier ° I 2 3 4 5 6 | 8 9 


Multiplicand into counters sa I 2 }a2] §$ | B31 S35 135 @ | Se 








If C,, denotes the value of the number standing in counter n, the sum 
required is 8C,+4C,+2C,+C, which may be obtained directly from the 
tabulator by suitable plugging. If the multipliers exceed the value 10, 
separate runs are required for each digit. 

The method is perhaps the easiest of all Hollerith methods to apply when 
only one sum of products is required and the multipliers are small, since no 
sorting is required. For a large number of related summations, however, 
it is uneconomical compared with other methods. 

(c) Binary self-rolling 

In this method the total standing in any counter at a change in the digit 
column effecting control is rolled into itself, thus doubling the amount 
standing in the counter. The method is particularly appropriate in work 


+ This is known in Hollerith usage as ‘the distribution method’. It is called binary 
distribution here to distinguish it from constant multiplier distribution. 
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such as Fourier analysis, where sets of factors have to be used with several 
functions to be analysed. As an example, suppose four functions a, ;, ¢;, d; 
are being subjected to Fourier analysis and that the first few values are 


t a b c d 
0 27 18 57 33 
l 35 12 51 45 
2 42 8) 40 47 


and suppose further that no value of any of the functions exceeds 63 (this 
limitation is necessary only for the purpose of this example). The set of 
cards for this work would be punched as follows: 





Field I Il III 
Binary Detail 
Code No. digit indicator Fourier 
Card No. f card |(one column only) coefficients 

I ro Yor 

2 02 Y X 

} O4 Yo 

5 O5 YXo 

of oI 

7 iI Yor Punch sine and cosine 
8 2 Yo factors relevant to 
) I X1 values of ¢ and the 
*) 14 XI trial periods being 
II 15 ° tested. 
12 Yor 
I I X1 
14 YI 
15 23 I 
10 24 YXor 
17 25 
18 2¢ Yor 





The code number of the card indicates the serial number of the factor 
(i.e. the value of the suffix ¢) and the binary digit in the factor. Card 15 is 
coded, 23 indicating that it applies to the data fort = 2and that the binary 
digit is the third lowest. The binary digit indicators, which are all punched 
in the same column, are only punched if the relevant binary digit is one and 
not zero; Y, X, 0, 1 refer to the digits in a, b, c, d respectively, e.g. 

a, 35 = 32+2+1 
or in binary form 100011, and hence cards coded 16, 12, 11 each have Y 
punched in them and cards coded 15, 14, 13 have not. The cards are sorted 
once and for all in descending order of the final column of the code number. 
In order to tabulate the summations involving the function a,, all details 
on cards punched with Y in Field II are entered into counters, using a dis- 
tributor to exclude details on cards not so punched. Control is on the final 
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column of Field I, counters being rolled into themselves at each change, 
The distributor is changed to X, 0, 1 for the b, c, d summations respectively, 
In this case, of data not exceeding 63, only six breaks in control would 
occur in each tabulation run. 

By using all the twelve lines of the cards, overpunching may be extended 
so that twelve functions appear on each set of cards. Although one card 
per binary digit is required, overpunching can effectively lead to very 
substantial card economy in suitable applications. The punching of binary 
digit indicators is readily done by reproduction from a library pack of 
binary coded cards. 

A variant of this method consists in punching all the binary digit indi- 
cators of each factor on one card in adjoining columns. This is advantageous 
when the details are being used with only a small set of functions. Its main 
virtue is that punching of the binary digit indicators is very much easier 
and faster than in the main method. In the long run, however, less card 
economy is achieved since Field II encroaches on the detail space and in 
addition, tabulation is slower since, during each tabulation, the plugboard 
has to be changed once for each binary digit in the functions, in order to 
obtain correct distribution. 


4. The National class ‘3000’ machine 

This machine is described with some applications by Comrie (4), Carter 
and Sadler (5), and Richards and Sadler (6). 

The process of rolling is performed by totalling, or subtotalling, the 
contents of a particular register with the carriage in the position for 
actuating the adding mechanism of another chosen register. Self-rolling 
requires two operations: firstly, the subtotalling of the contents of a register 
with the carriage in that register position, and secondly, the entering of the 
printed subtotals into the same register. Distribution may be achieved by 
the operator choosing the registers into which a number has to be entered; 
on occasion this can be done by using multiple stops. 

All the methods of Hollerith summation, except binary distribution, 
involve sorting data into certain orders. Some means of doing this is 
required before the methods are applied on the National machine. A very 
effective way of sorting data involves the use of Paramount cards. These 
cards, which are provided in various sizes, are punched with a series of holes 
along the edges. Each hole represents a certain number, letter, or value, 
according to the code adopted. After the necessary data are entered on the 
face of the card, a pair of ticket nippers can be used to slot out the holes 
which represent the particular details governing the orders in which the 
cards will be required. When a group of cards has been slotted, it is simple 
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to arrange them in any prescribed order by passing a needle through the 
appropriate hole and lifting; the required cards then drop out, being thus 
separated from those remaining on the needle. 
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Fic. 1. Paramount card prepared for use in triple product summation on National 
Class ‘31’ accounting machine. 

For computing purposes, binary coding of the numbers slotted is best 
for two reasons: firstly, the numbers are required in binary form in two of 
the methods to be used, and secondly, the binary code is obviously the 
most efficient coding for this system of sorting, which depends on two 
‘symbols’ only, slotted and unslotted holes. 

A Paramount card, coded for use in a triple product summation, is 
illustrated in Fig. 1. 

For many purposes, however, only one initial sort is required, and in 
such applications any simple type of filing card, or even slip of paper, may 
be adequate, the sorting being done by hand. 

The conventional abbreviations used by Carter and Sadler to describe 
machine operations are used here. The symbol C,, is introduced, however, 
to indicate the number standing in register n, . the operation 7', or ST’, 
will result in the printing of the value ¢ 

(a) Binary SE lf-rolling 

Suppose the sums of products 


. : 
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are required. On the body of the ¢th card is written 


11 
Mats Agty--- Xa ays 2% 
v 


At a convenient coding position wu, is slotted. If more than eleven sums are 
required, the relevant «’s may be written on the same set of cards. If more 
than one function wu, is being used, as may occur, for example, in Fourier 
synthesis or analysis, then the different u,’s may be slotted until coding 
capacity is exhausted. 

The machine set-up and operation sequences are shown in Table 1. The 
cards are sorted on the column of the highest binary digit of u, and those 
slotted in this column are tabulated. After the self-rolling sequence the 
cards are re-sorted and those slotted on the next highest binary digit column 
are tabulated. This process is repeated until the end of the tabulation of 
digit 1 is reached, at which the final sequence is to total all registers giving 
the first six sums of products. Repetition gives the next five and }-check. 
If this check does not give agreement the error can be located by checking 
subtotals at self-rolling sequences and finally by checking the printed 
record. The result can be adjusted without complete reworking. 


(b) Successive digiting 

The set-up for this method, which is directly analogous to the Hollerith 
method, is obvious, subtotals being formed in three registers and rolled into 
the other three. Tabulation may proceed from highest to lowest multiplier 
in one run, but if many multipliers are missing it may save time to split 
the higher multipliers and tabulate their associated multiplicands with 
otherwise blank multipliers (e.g. 68 might be tabulated as 50 and 18), filling 
all blanks below a chosen maximum. In this method errors are more easily 
detected than in the binary self-rolling method, but only half as many sums 
are formed at each run. Moreover, there are more control sequences. In 
practice the choice between these methods in a particular case depends on 
the number of sequences required: this can be readily estimated before the 
start of tabulation. 


(c) Manual binary distribution 

Like the identical Hollerith technique, this method is not practicable if 
multipliers are large, unless totals are small enough to permit of more than 
one being accumulated in each register. 


(d) Constant multiplier distribution 
If, as in numerical quadrature by Cotes-type formulae, there are a few 
constant multipliers, distribution may be effected by allotting a register 
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to a specific multiplier. For example, the repeated application of the 
six-strip integration formula (5) of section 2 can be achieved by the scheme 
Stop +8 44 48 464,245 
Sequence of tabulation to Ss Se Js Se Ss 
te Ji Ss 


etc. 
Then 


7 Ww Y sy ons ore Y 
JA) - 140 (82¢ 1 +216€ 37 270; T 272C,—4(fo4 Son)}- 


« 
0 


This is of course identical with the desk machine method, but the printed 
record provides a visual check, whereas the desk machine method can only 
be checked by repetition. If common factors can be found, the factors need 
not be limited to six, e.g. 8, 9, 12, 15, 17, 23, 29 might be distributed with 
seven stops +1, +2, +3, +4, +1+2, +144, +1+42-+3, and the sum- 
mation would be 8C,+-9C,+-12C,+ 15C,. 


5. Negative numbers 

In Hollerith work negative numbers may be entered into counters in 
complementary form. The three-digit numbers — 144, +144 are punched 
in four columns of the card as 9855 and 0144 respectively. Multiple plugging 
of the thousands digit to all higher digits of the counter ensures that the 
correct number enters the counter. The 9’s complement 9855 is used 
instead of the true complement 9856 because of the ‘elusive one’ feature 
of the tabulator. If preferred, this operation may be avoided either by 
immobilizing the elusive one relay, or, provided totals are not expected to 
exceed say seven or eight digits, the 9 indicating the complement need not 
be plugged through the whole counter; in either case the true complement 
isused. Whenever possible, it is simpler to make all data positive by working 
from a false origin. For single summations all multiplicands can be punched 
positively, and the sign of the multiplier changed accordingly and punched 
is X for positive, Y for negative. It is then possible to sort and tabulate 
the positive and negative products separately. Alternatively, if counters 
are available, sorting is avoidable by distributing multiplicands into 
separate counters ac -ording to the sign of the multipliers. 

Negative numbers are not so troublesome on the National machine. 
Except when registers are doubled up, the multiplicand can be added 
positively or complementarily into the registers according to the sign of 
the multiplier, the operator doing the conversion to complementary form 
is necessary at sight. When Paramount cards are used, multipliers may 
be written on the cards as well as slotted, using red and black figures for 
negative and positive numbers. This gives a good visual aid to the operator. 


5092.26 
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6. The National class ‘31’ machine 

This new machine differs from the class ‘3000’ in two important respects 

it has ten registers (numbered 1, 2,..., 9, 2) and it is possible to enter 
negative numbers directly into all ten registers. This increase in capacity 
and flexibility is enhanced by the further feature that, by the use of 
additional control stops, complete sequences of control operations can be 
carried out fully automatically. The machine can pass over control stops 
on tabulating sequences and vice versa. The increased automaticity not 
only greatly increases the speed but also the accuracy of operation by 
reducing the possibility of operating errors. 

The following two examples, which have been tested, demonstrate the 
way in which the greater flexibility and capacity can be used to improve 
on class ‘3000’ methods. 


(a) Sums of products by a self-checking binary self-rolling method 

The set-up is givenin Table 2, cards being prepared as in section 4 (a). 
On tabulating, sequences «)...1g are entered into registers 1...9 respectively 
and printed. At the same time each enters register x. The negative entry 


9 
of >} a, on the —x stop thus reduces the contents of register x to zero 
i=1 


and the “Total 2’ stop then clears register x and prints the total C,, = 0. If 
an error has been made this is shown immediately as a non-zero entry in the 
column of print. 

At changes in the binary digit, on the first control sequence, the stop 
ST.,,, subtotals register 1, prints C, and enters C, into register x; the 


+X 


grees Ug 


9 
stop ST’, then gives 2% which is printed. Thus, after entering C),...,C 


i= 


9 
on the second control sequence the sum standing in register x is > C; and 
i=1 


the operation set C, then clears register x and the ordinary tabulating 
sequences may be proceeded with. 

By having this automatic visual line-by-line check, all errors, which are 
usually errors of setting, can be corrected as they occur. 


(b) Summation of triple products 

The summation of ¥ xyz can be performed by combining binary self- 
rolling and any distribution method. Consider the sum in the case where 
x;, Y;, and z; take the values 


D ] 2 3 t 5 
x; — 39 24 53 19 31 
Yi 2074 562 1093 —1762 —1891 


—523 —898 262 
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Here we see that x, requires six binary digits, and z; ten. Assuming that 
the rest of the series require the same numbers of digits, this could be 
performed on the ‘3000’, distributing on 2; and self-rolling on z,;, but there 
igmuch repetition in the tabulation. On the ‘31’ it is more easily performed 
by distributing on z, and self-rolling on 2;. It is first necessary to make 2; 
ind z, positive by transferring signs to y;. Cards would be prepared as 


follows (see Fig. 1): 


Card number | 2 3 4 5 

r. (slotted) 1.2.3.6 4,5 1, 3, 5,6 1, 2,5 1, 2,3, 4,5 

j, (written on 2074 562 1093 1762 — 189] 
card L3.4,2 2&2 2,3,9 1,4,5,6,7 2,3,4,6,9 


The cards are sorted on the highest 2 digit, tabulated by binary distribution 
of y in the appropriate registers according to the coding of z; and the 
registers are all self-rolled at each digit change in x. The tabulation of the 


summation given here starts as follows: 


Binary 
ligit Card 
1 No. C; C. C, c 
6 l 2074 2074 ... 2074 
3 1093s... — 1093 
2074 | ar 1093 2074 
2074 981 an 1093 2074 
5 2 562... 562 — 562 
3 1093s... — 1093 
4 1762 soe CHC. 


An application of binary self-rolling combined with constant multiplier 
distribution is to the evaluation of integrals of the type [ f(x)g(x) dx where 
1 Cotes-type quadrature formula is applicable, particularly when f(z) and 
(x) can be written on Paramount cards directly from tables. 


7. Conclusion 

Direct comparison of the relative economies of the methods described is 
difficult since such a comparison is bound to depend largely on the cireum- 
stances of the individual. Information on the relative times required for 
various methods is given in Appendix II. 

A comparison between the National and desk machine methods can be 
made most easily by considering the case of the progressive digiting method. 
Provided that there are few missing multipliers, a single sum of products 
can be formed, using two registers on the National Class ‘3000’, at much 
the same speed as the multiplicands can be added on an electric desk 
machine. Thus the initial times taken on the National ‘3000’ and on electric 
desk machines to form a single sum of products is about the same. However, 
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on the desk machine method there are three main sources of error, in setting 
the multiplicand, in multiplying, and in recording intermediate totals, 
whereas only the first occurs in the National method. Moreover, the incor- 
porated checks enable mistakes to be located and corrected on the printed 
National record without reworking. On desk machines, although }-checks 
show up the occurrence of mistakes, they can only be located by reworking 
all the summations from the stage at which the check of the intermediate 
record shows the working to be correct. This time saving, even with only 
one error in several summations, is considerable. 

As far as cost is concerned, a rough rule (based on charges made by the 
respective makers of the equipment) is to take the cost of Hollerith tabulat- 
ing time at four times that of National machine time, and then to add the 
cost of cards and their preparation—including the clerical cost—and also 
the cost of any post-mechanical work such as desk combination of results, 

From extensive work performed by the methods described in this paper, 
the author is of the opinion that the National ‘3000’ methods are slightly 
more economical on this basis than the comparable Hollerith methods for 
computations of moderate extent, i.e. computations from the stage at 
which desk machine methods need to be superseded until the stage at which 
time saving is the dominant factor. At this stage the Hollerith methods 
supersede the National ones. The National methods do, however, make 
more exacting demands on the operator’s concentration. 
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APPENDIX I 


Hollerith Card Preparation for Serial Correlation Evaluation by 
Progressive Digiting Method 


‘or three-figure data: columns 1—39 of cards are Control Fields 1-13 each of 3 
lumns; columns 40—72 are Detail Fields 14—24 each of 3 columns; columns 73—76 
ive the ¥-check Field 25 of 4 columns. The remaining columns are used for the 


rial number of the card. 


Cards are punched as follows: 


ld | ms . 2 “a i6...23 24 25 No. 
urters Uy 1] 
later uy Us x10 
lisearded) Uh, Ue... te Uy, xX] 
12 
U, Me Me __ Slre > u; l 
7 
23 
? A 
{ ba Msg tag. - thes Uos u; | 12 
13 
144 
PY 
é( Ujo + + + Uyo2 Uy33 Uiza Uig5-++Ui4a3 0 Ua at 133 
134 
11 
Un—132 “n_-121 Ln Unsi Un+e Uni Ujsn n 





The punching is done by (i) hand punching Fields 24, 25 and serial number on each 
ard; (ii) off-set gang punching Fields 13—23 starting from 24; (iii) ‘off-set reproducing’ 


fields 1-12, i.e. reproducing fields 24 into 12, 13 into 11, 12 into 10... 3 into 1 from 


ard number n to card number -+- 22. This may be done in runs of 22 cards; on com- 


pleting each run the 22 newly-punched cards are transferred to the feed to reproduce 


» 


the next 22, and so on. Without the preliminary ‘starter’ cards, the transferring 


would be necessary after every run of 11 cards. 


It is convenient to have the detail fields together with the }-check field, as a 


multiple of 3 in rolling tabulation or 6 in the summary punching method. The set-up 


lescribed gives the maximum number (143) of sums of products obtainable on the 
ard subject to this condition. For two- or four-digit fields, divide the available fields 
nto Controls (C’) and Details (D) so that CD maximum, subject to C+D = con- 


stant, and so that D+ 1 is a multiple of 3 or 6. The cards are, of course, tabulated 


n control fields successively. The }-check reconciles successive runs of 11 sums of 


products. 


This layout gives many more sums of products than will usually be required. It is 
ore economical than Hartley’s original layout whenever the original requires 


reproduction of two or more pack of cards, since in that case the new layout requires 


fewer cards and, at worst, no more reproducer runs than the original. 
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APPENDIX II 


Comparison of Times taken by Different Methods 


There are two types of problem, in one of which the data is already tabulated and 
card preparation can be started immediately. In this case the preparation takes about 
the same time for either Hollerith or Paramount cards. It will often be possible to 
list the data on the National machine at the same time as the sums are being formed 
for >-checks. The lists, cut up appropriately, may then be pasted on to the sorting 
cards. In the other case, where data has to be collected from several sources, as, e.g, 
statistical data from several yearly reports, Paramount cards may be quicker to 
prepare than Hollerith cards, since the cards themselves may be used as work sheets, 
and when the data is in the final form ready for use, the Paramount card is re» dy, 
whereas the Hollerith card has yet to be punched. This affects the decision a: to 
which machine to use, since the time saved in card preparation may offset that lost 
in tabulating. 


(a) Serial correlations of short time series 

Sums of products for 11 serial correlations and }-checks of 72 consecutive monthly 
values of the products of inertia of the Earth calculated from atmospheric pressure 
were evaluated. No clerical preparation of data was required, and card preparation 
took about the same time for both types of card—for 72 cards it was not considered 
worth while setting up special plugboards for the reproducer in the Hollerith method. 
The times taken were: 


Hollerith method for forming the 11 sums and check from 72 cards: 


Progressive digiting and rolling . ; . 40 minutes 
(Subsequent desk time for combination . 10 minutes.) 


National ‘3000’ methods for the same work: 


Progressive digiting and rolling . : . 140 minutes 
Binary self-rolling  . ; ' ; . 125 minutes. 
(No desk combination required.) 


It will be seen that, equating the cost of 4 hours of Hollerith to 1 hour of National 
tabulating, the National processes are more economical in this case. It appears from 
other computations carried out that for sums of up to about 200 products the National 
process is economical, but beyond that the Hollerith is to be preferred for its lesser 
demands on the operator and for the time saved. 

No account has been taken of the times taken to set up the machines, as in long 
calculations these ‘overheads’ become insignificant. 


(b) Periodogram analysis 

The extraction and editing for the first time of sine and cosine factors in the permuted 
forms required for periodogram analysis is a fair example of a calculation in which 
clerical preparation of data should be estimated. 

With the same two series as above, 54 sets of factors for 27 trial periods were taken 
from tables and the sums formed for }-checks. The extraction was done straight on 
to 72 Paramount cards in about 4 hours. The Hollerith ecards (180 in number) took 
a further 1} hours to prepare. 
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For this work (54 sums of products+5 checks, average length, 60 products), the 


tabulating times were 


Hollerith: 
Binary self-rolling . : ; . . 40 minutes 
(Time for preparing extra binary cards about 90 minutes.) 
Progressive digiting and rolling . : . 75 minutes 
(Desk machine combination ‘ ; . 45 minutes). 


National ‘3000’: 
Binary self-rolling . ‘ . . . about 5} hours. 


It is seen that in this case there is little to choose between the methods apart from 
the initial advantage the National machine methods have for quicker card prepara- 
tion. It is noticeable, however, that if many Fourier analyses of syntheses are made, 
the one set of Hollerith cards with the factors punched can be used over and over 
again, using overpunched binary coding for the function being analysed or synthesized. 
If this is done by reproducing from a binary pack, the binary self-rolling method 
becomes not only the quickest in actual time but the most economical of all from the 


urd-usage point of view as well. 











ON THE TRANSFORMATIONS OF SINGULARITIES 
AND LIMIT CYCLES OF THE VARIATIONAL 
EQUATIONS OF VAN DER POL 


By A. W. GILLIES (Northampton Polytechnic, London) 
[Received 6 August 1952] 


SUMMARY 

The investigation of the solution of van der Pol’s equation with forcing term leads 

to equations for the amplitude, 6, and phase of the oscillation, ¢, of the form 
b b(1—b*?)— F cosd, 
bd —ba+ F sind. 

The solution of this autonomous system of first-order equations has been dis- 
cussed by Cartwright. By considering the isoclines on the plane with (6,4) as 
polar coordinates, it is shown that Cartwright’s solution is incorrect in one range of 
the parameters. The corrected solution is given. In consequence of this, it is shown 
that the hysteresis effects to be expected for a van der Pol oscillator with increasing 
and decreasing frequency are confined to a narrower frequency interval and are less 
varied in character than was suggested by Cartwright’s solution. 


Introduction 
THE investigation of the solution of van der Pol’s equation with forcing 
term (1), (2), (3) leads to the study of the variational equations of the form+ 


b = b(1—b?)— F cos d, (1) 

bé = —ba+Fsind (2) 

in which 6 and ¢ denote respectively the slowly varying amplitude and 
phase of the forced oscillation (dots denote time derivatives). Equations 
(1) and (2), in which } and ¢ are polar coordinates (so that the vector (b, 4) 
is the representative vector of the oscillation) form an autonomous system 
of first-order equations, the singular points of which represent the harmonic 
solutions of the original van der Pol equation, and the limit cycles of which 
represent combination oscillations of the van der Pol equation. The 
present paper is concerned with the variation in the number and positions 
of the singularities and limit cycles of this sytem when the parameters F 
(representing amplitude of applied e.m.f.) and x (representing the frequency 
of applied e.m.f. measured from the free oscillation frequency) are varied. 


+ With the equation 

d2yv vy dv 

en Ss adele aioe ee Sot 

i €wo(1 73) = wot Ba?* sin wt 
and w = wy+Aaw, write v = 2v, b cos(wt+¢), equations (1) and (2) are obtained by assuming 
b and ¢ to vary slowly, and setting F = (Bw)/(2ew, v9), c = (2Aw)/ (ews), and 2t = ewy7- 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 2 (1954)] 
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The number and distribution of the singularities is already well known. 
it is usual to represent them by plotting y = 6? for the singular points 
against x for different values of F which give resonance curves for the 
van der Pol oscillator having the equation F? = y|(1—y)*+2?]. These are 
shown in Fig. 1. For |a 1/\3 there may be one, two, or three singular 
points; for a! > 1/\3 there is one. The type of singular point is also 











—1 0 +1 x 
Fic. 1. van der Pol resonance curves. 


ee l 
Sve ge" Ge Bee Ben 26 


The curves correspond to F? = 

The broken lines separate regions for which the points 

represent different types of singularity as indicated. The 

portion of the figure within the rectangle is reproduced 
on an enlarged scale in Fig. 7. 


indicated on the figure; those represented by points within the ellipse é 
having the equation (1—y)(1—3y)+a? = 0 are cols; those represented by 
points outside & and above y = 3 are stable nodes or foci; those outside 
&€ and below y = 4 are unstable nodes or foci; points above y = |x 
represent nodes, those below y x| foci. 

The points on & represent second-order singular points which separate 
into two simple singular points on one side and become imaginary on the 
other side. The extremities of the axis of & parallel to the x-axis, |x| = 1/v3, 
y = 2/3, represent third-order singular points when all three singular points 
are coincident. 
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It is also well established that for F* < 4/27 there is no limit cycle so 
long as three singularities exist, the solution curves running from infinity 
or from the unstable point to the stable point, apart from the two which 
run into the col.+ When the stable node and the col move into coincidence 
and disappear, at the point where the resonance curve crosses &, a limit 
cycle forms from the separatrix running from the col to the stable node, 
and we are left with a stable limit cycle surrounding an unstable node or 
focus. 

Furthermore, for any F and |x| sufficiently large there is a single limit 
cycle surrounding an unstable focus, and for F? > 8/27 there is never 
more than one singular point; when the singular point crosses from 
y <4(orb < 1/v2) toy > 3 (orb > 1/2) it changes from an unstable to 
a stable point by the stable limit cycle shrinking up to the singular point 
and disappearing. In fact it may be proved by the standard procedure (4) 
that over the whole length of y = } outside @, i.e. \v| > 4, a singularity 
crossing y = 4 for |x| decreasing becomes stable by a stable limit cycle 
shrinking up to it and disappearing. 

The question at issue is the transition for F? < 8/27 from a stable limit 
cycle surrounding an unstable focus for |2| large, to no limit cycles and 
either a single stable node, or a stable node, a col, and an unstable node or 
focus, for || small. 


Representation on an (x, F) plane 

We are concerned with the pattern of integral curves on the van der Pol 
plane (on which 6 and ¢ are polar coordinates) for any particular pair of 
values of F and x, and with the way in which this pattern changes when 
the parameters F and 2 are varied. For this purpose the resonance curves 
on the (x, y)-plane are not convenient, and a representation on a plane on 
which x and F are taken as coordinates will be found more useful. On this 
plane each point (x, F) represents a particular pattern of integral curves on 
the (6, ¢)-plane. Continuity requires that if x and F are varied continuously 
the pattern of integral curves must also vary continuously. Thus, suppose 
a particular point (x, F) corresponds to a pattern having, say, a single 
stable limit cycle surrounding an unstable focus, then we can draw a circle 
centre (x, F) and suitable radius p such that every point (x-+-8x, F-+5F) 
within this circle also gives a pattern having one stable limit cycle sur- 
rounding an unstable focus. Thus particular types of integral curve pattern 
will correspond to points within certain regions of the (2, F)-plane, these 
regions being bounded by curves, the points of which represent transitional 


+ See Appendix. 
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patterns. Thus, the region containing points (x, F) corresponding to pat- 
terns such as those just mentioned, havingone stable limit cycle surrounding 
an unstable focus, will have a part of its boundary formed by a curve whose 
points represent transitional patterns for which the limit cycle has shrunk 
up toa point. Points on the other side of this boundary will correspond to 
patterns having no limit cycle and one 
stable focus. The radius p will tend to 
zero as the point (x, F) approaches such 
a boundary. 

Clearly there are many different ways 
in which the pattern of integral curves 
may be characterized with corresponding 
subdivisions of the plane into regions 
bounded by curves representing as many 
different transitions from one type to 
another. Fig. 24 shows some of these 
curves embodying the results already 
summarized in the introduction. (The 
whole figure is symmetrical about the 
F-axis and only the half for which z is 
positive is shown.) 

The region OA B contains points repre- 
senting patterns having three singular 
points, points on the boundaries repre- 
sent patterns containing second-order 


singular points. On OA, given by 


F? — (2/27)[1+ 92?—(1—3z2?)3?], 








the double singular point corresponds to 





the upper boundary of the ellipse &, while 
Fic. 2. Diagrammatic represen- 


m AB, given by . 
tation of the (x, F’)-plane. 


F? (2 27)[1 1 Gy +-(] 3a?) 3/2) 
the double singular point corresponds to the lower boundary of &, these 
equations being obtained by eliminating y between the equation of the 
resonance curve and that of &. The cusp at A (1/¥V3, ,/(8/27)) represents 
the pattern containing the triple singular point at x = 1/v3, y = 2/3 on 
the figure showing the resonance curve. H is the point on OA with 
ordinate ,/(4/27). 


The hyperbola F? = }(x?+-}), shown as DCK M, contains points repre- 
senting patterns having a singularity on y = orb = 1/v2. It touches BA 


at the point C(4,4), which represents the pattern having a second-order 
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singularity at y = }, i.e. it is represented by the intersection of & with 
y = 3 0n the resonance curve figure; K is the point on this hyperbola with 
ordinate ,/(8/27). 

The point L (v5/4,3/(4v2)) at which the hyperbola intersects OA, on the 
other hand, represents a pattern having a double singularity corresponding 
to a point on the upper boundary of & with the third singularity on y = 1. 

For points on the are DC of the hyperbola the singularity on y = } is 
a col within the ellipse &. For points on the are CLK M the singularity on 
y = } is a centre transitional between a stable and an unstable focus. On 
the right of this arc the focus is unstable and on the left stable, a stable 
limit cycle having shrunk up to a point and disappeared as the arc is 
crossed from right to left, with decreasing x and fixed F. 

Ignoring the are DC of the hyperbola the plane is divided into four 
regions. In region I above the line BCALKM there is a single stable 
singularity which is a node or focus, and there are no limit cycles. In 
region II to the right of OL KM there is a single unstable focus surrounded 
by a stable limit cycle. In region III, within OLCB, there are three 
singularities, a stable node or focus, a col, and an unstable node or focus. 
In region [V—the small curvilinear triangle CLA—there are three singu- 
larities, of which two are stable and the third is a col. 

Clearly one could add to the figure by adding the locus on which a 
singularity is transitional between a node and a focus, or in various other 
ways, the plane being broken up into smaller regions, in each of which the 
type of integral curve pattern is more completely characterized. This is, 
however, unnecessary for the present discussion. 

To each point on the (x, F’)-plane corresponds a pattern of integral curves 
of equations (1) and (2) with parameter values equal to the coordinates of 
the point. It will save circumlocution if the point is identified with the 
corresponding integral curve pattern so that a statement such as ‘a limit 
cycle through the col occurs at a point P’ is to be understood to mean that 
‘a limit cycle through the col occurs in the pattern of integral curves which 
corresponds to P’. 

If we follow a resonance curve on the (x, y)-plane we vary x and keep 
F constant, and thus follow a parallel to the x-axis on the (x, F)-plane. 
Moving from large values of 2 towards 2 = 0, i.e. from right to left 
on the (x, F)-plane, we have, for F > ,/(8/27), on the extreme right, 
a stable limit cycle surrounding an unstable focus. Moving to the left the 
limit cycle shrinks up round the focus, disappearing when we reach the are 
KM, the focus then becoming stable. Continuing the movement to the left, 
the focus changes to a node on a locus not drawn, before we approach the 
maximum of the resonance curve on the F-axis (x = 0). 
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For F < ,/(4/27) we have again on the extreme right a stable limit cycle 
surrounding an unstable focus, this being in fact the situation on the 
extreme right for any F. This time, when we reach the are OH of the curve 
0A, a second-order singularity forms on the limit cycle. To the left of OH 
the second-order singularity separates into a stable node and a col, the 
limit cycle ceasing to exist, and giving place to a separatrix joining the col 
tothe stable node, which separates integral curves from the unstable focus to 
the stable node from integral curves coming in from infinity to the stable 
node. 

For the region ,/(4/27) < F < ,/(8/27) we still have a stable limit cycle 
on the extreme right and no limit cycle near to the F-axis on the left, so 
that somewhere in between the limit cycle must disappear, and during this 
process the requirement of continuity must be satisfied both above, where 
the limit cycle shrinks up to a point, and below, where the limit cycle is 
broken by the appearance of the second-order singularity. 

We cannot continue the two transition loci MK and OH to meet at L, 
since at L we cannot have a limit cycle which has shrunk up to a point 
round one singularity, and at the same time passes through a quite separate 
second-order singular point. The problem then is how the transition is to 
be effected in this range while satisfying the requirements of continuity 
throughout the range. 

Two possible ways will be described in which the transition can be 
accomplished which satisfy the requirement of continuity and which con- 
form to the other facts regarding the differential equations (1) and (2) which 
have been mentioned. The first of these is the solution proposed by Cart- 


wright (3), (5). 


Cartwright’s solution 

The critical region of the (x, F’)-plane is the region containing the points 
CALKH. This region is drawn in Fig. 24 on a distorted scale so that the 
different regions may be shown more clearly. The information obtained so 
far is represented diagrammatically by representing the limit cycles by 
horizontal lines. In the upper part of the figure these terminate on MK 
and in the lower part they terminate on OH. In between the precise way 
in which they terminate has yet to be decided. 

Still considering the sequence of integral curve patterns for fixed F and 
decreasing x, i.e. moving from right to left on a line parallel to the x-axis on 
the (x, F)-plane, Cartwright follows the way in which this sequence may 
change, having regard to continuity, when the line is moved downwards 
starting with F large. For a fixed F > ,/(8/27) the limit cycle shrinks up 
toa point when the transition locus M K is reached as x decreases. Therefore, 
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Cartwright says, this must continue to be the case right down to F = 4} when 
the line parallel to the x-axis passes through C. For a fixed F between } 
and ,(8/27) the limit cycle shrinks up to a point when CK is reached. 

Thus a limit cycle shrinks up to a point at C for F = }. For a fixed F 
slightly smaller than } and to the right of C we must have a limit cycle 
which is shrinking as x decreases, approaching the value x = 4 of C; while 
for the same F and to the left of C no limit cycle can exist. Therefore the 
limit cycle must disappear by passing through the col at a point near to ( 
within region III. This is possible with a limit cycle which is approaching 
vanishingly small because the col and unstable focus are near to each other 
approaching to coincidence at C. 

We are thus led to a new transition locus CN representing patterns in 
which the limit cycle passes through the col. This locus must extend from 
C to meet OHA at the point N above H at a value F = ,/(4/27+-8), where 
6 is not determined. 

For a fixed F below this value the limit cycle disappears as x decreases, 
by passing through the second-order singularity in the pattern represented 
by a point on HN, and becoming a separatrix. The requirements of con- 
tinuity are thereby satisfied at both extremities of the range. This solution 
is shown diagrammatically in Fig. 2B. 


The alternative solution 

We again consider the sequence for decreasing x and a fixed F, i.e. 
following parallels to the z-axis from right to left, but this time for succes- 
sively increasing values of F. Commencing with F < ,/(4/27), we have a 
limit cycle which ceases to exist when OH is reached and the second-order 
singularity appears. 

For each fixed F the limit cycle may continue to disappear in this way 
until we reach the parallel to the x-axis through C, i.e. F = 4. For a fixed 
F > 4 the focus changes from an unstable to a stable one as x decreases 
by a stable limit cycle, shrinking up to it on CLK; and if the original limit 
cycle has already disappeared at a greater x on SL (S being the point on 
OHA with ordinate }) by passing through the second-order singularity, 
another must arise, and it can only do so by forming through the col and 
then shrinking up round the focus. This can occur without violation of 
continuity again because the col and unstable focus coincide in the pattern 
represented by C. 

We are thus led to postulate another transition locus CP representing 
patterns for which a limit cycle forms through the col as x decreases. If this 
locus meets OHA at P, then for a fixed F = Fp and decreasing 2, the limit 
cycle passes through the second-order singularity at P, but instead of 
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degenerating to a separatrix when the singular points separate, it continues 
to exist, shrinking up round the unstable focus to vanish on CLK, when the 
focus becomes stable. 

For a fixed F > F, and decreasing x, the second-order singularity 
:ppears outside the limit cycle, the latter continuing to exist until it shrinks 
up to a point on CLK. P must be below L, since for a fixed F > F, and 
decreasing x the focus becomes stable before the col appears. 

Thus we have an alternative sequence of transitions which satisfies the 
requirements of continuity. This is shown in Fig. 2c. 

Comparison of the two solutions 

The two proposed solutions give the simplest ways in which the transition 
could be accomplished. The possibility of additional limit cycles appearing, 
or of the original splitting and recombining, has not been formally shown 
to be impossible, but even if this happened we must still have a limit cycle 
forming or disappearing by passing through the col in some such way as has 
been indicated. 

When considering a system of equations dependent on a single parameter, 
the formation of a second-order singularity on a limit cycle is normally a 
prelude to the degeneration of the limit cycle to a separatrix when the 
parameter is further varied. The persistence of the limit cycle after the col 
and node have separated would be exceptional. It may have been such a 
consideration, or merely the fact that for a fixed F fewer transitions are 
required, that led Cartwright to the first of the proposed solutions rather 
than to the second. 

When we have a system dependent on two parameters, however, the 
persistence, or not, of the limit cycle is dependent on the way in which 
the parameters are chosen. With Cartwright’s solution in Fig. 2B, at 
the point V the second-order singularity has formed on the limit cycle. 
If new parameters €, 7 are chosen so that » remains constant on QN R, 
then for variation of € the limit cycle which exists on QN passes through 
the second-order singularity at NV, but continues to exist when the second- 
order singularity has separated into node and col on NV R. 

Thus it would appear that there is no essential difference as regards 
simplicity between the two proposals. In any case, simplicity does not 
prove correctness and it is a matter for further investigation whether the 
solutions of the system (1), (2) behave in either of these ways or in some 
more complicated manner. 

Discrimination between the two possibilities 


The system of equations (1), (2) cannot be integrated by known methods, 
and there are no general theorems regarding the formation of a limit cycle 
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through a col which would enable a decision to be made between the 
different possibilities which have been described. We are therefore thrown 
back on numerical methods, which can prove very laborious, and possibly 
unreliable in critical regions. 

However, it has been found that an examination of the isoclines for a 
comparatively small number of particular cases enables a decision to he 
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m, Mz 
Fic. 3. Integral curve confined between 


two polygons. 


reached. Referring to Fig. 3, if we have a number of isoclines corresponding 
to gradients m,, m5, Mg,..., such that nowhere in this region does an isocline 
have a point at which the gradient of the isocline is equal to the corre- 
sponding m (so that the integral curves have no point of inflexion in this 
region), then an integral curve starting from a point A on the m, isocline 
must lie between two polygons AB, C, D,... and A B,C, D,..., the vertices 
lying on the successive isoclines, and such that the gradient of each side of 
the first polygon is that of the isocline through its vertex nearer to A, while 
the gradient of each side of the second polygon is that of the isocline through 
its vertex farther from A. This is shown in Fig. 3. 

A critical case is that corresponding to the point S (Fig. 2) having the 
same F as C, i.e. F = 3, and the value of x which gives the double singu- 
larity, namely 2 = 0-5221. A number of isoclines for this case are con- 
structed in Fig. 4, the gradient being marked on each curve. The curves 
6b = Oand ¢ = 0, marked T and R respectively, indicating transverse and 
radial directions respectively, are also shown. All these curves pass through 
the singular points, the R and 7 curves being tangential at the double 
singularity and the infinite gradient curve having a double point at the 
same point. 

Reference to Fig. 2 shows that, if uhe first solution is correct, a limit cycle 
enclosing only one singular point should have expanded with increasing 2 
from a point at C; and at S the double singular point should be outside this 
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limit cycle. If, on the other hand, the second solution is correct, a limit 
eyele through the double singularity should just have formed at S. In the 
first case the integral curve leaving the double singularity will wrap itself 


round the stable limit cycle. All integral curves finishing at the double 














Oo 


05 b, 
Fic. 4. Isoclines for F = 4,x = 0-5221 on the van 
der Pol plane, polar coordinates (b,¢), or Cartesian 
coordinates (b,,b,). Only one quadrant of the 
plane is shown; the heavy lines represent integral 
curves. 


singularity do so along the common tangent to the 6 = 0 and ¢ = 0 curves, 
except for two which approach the singular point parallel to the b, (ordinate) 
ixis.} These two would be separatrixes separating curves approaching the 
limit cycle from those ending at the singular point, as indicated in Fig. 5 (vii). 

In the second case, the curve leaving the double singular point must 
return to it again as shown in Fig. 5(iv). The isoclines of Fig. 4 show clearly 


+ Equation (3) of the appendix gives the directions through the singular point on putting 


R= 0and ® = 0, giving sin 20 = x/y. If one direction is parallel to the b,-axis @ = 47—¢p, 
€. sin 24, x/y), and on substituting sin dy = by 2/F, cos do bo(1—b3)/F, this reduces 
to (l—y,)(1—3y,) +2? = 0. 
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Fic. 5. Sequence of changes in integral curve pattern. 
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The critical region is between F = 3} and F = ,/(9/32). For 
increasing x we have first of all the single stable node S. The 
col C and stable node s form on the lower boundary of 4, 
the node s changing very quickly to a focus to give Fig. (i). 
The focus moves inwards towards the origin and becomes 
unstable (U), throwing off a stable limit cycle giving (ii). For 
} < F < & the limit cycle expands to pass through the 
col (iii). The curve leaving the col on the right takes a pro- 
gressively wider sweep, passing to the other side of S before 
S and C move into coincidence in (iv). The double singularity 
then disappears leaving a new limit cycle surrounding the 
unstable focus U, (v). 
For F = Fp the expanding curve leaving the col just reaches 
C again when C and S become coincident, running into C in 
the direction parallel to the b,-axis, the sequence being (i), 
(ii) to (vi), the second-order singularity disappearing to give 
(v) again. 
For Fp < F < 4/(9/32), C and S unite before the limit cycle 
has expanded sufficiently to reach the col, the sequence being 
(i), (ii), (vii), to (viii), the disappearance of C and S no longer 
producing a new limit cycle. 
For F > 4/(9/32), C and S unite and disappear before s has 
become unstable and no limit cycle forms until s becomes U. 
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that the second of these is correct. The integral curve leaving the double 
singular point is on the right of the polygon abcdefg, the first side ab extend- 
ing from 7' to the —4 isocline with gradient — 4 and tangential to T ata, and 
each side having the gradient of the isocline on which it ends. The integral 
curve entering the double singularity parallel to the b,-axis from below lies 
to the right and above the polygon stuv each side of which has the gradient 











Fic. 6. Integral curve pattern for F? = 6-625/27 and 
x 0-5166. 
The value of F® is slightly less than } and 2 such that 
the double singularity corresponding to the intersection of the 
resonance curve with the upper boundary of the ellipse & occurs. 
The limit cycle has formed through the double singular point. 
In addition to the integral curves the ¢ = 0 and 6 = 0 curves 
are shown marked R and T respectively. 


of the isocline on which it starts, the first side being parallel to the b,-axis. 
These are clearly consistent with Fig. 5 (iv) but inconsistent with Fig. 5 (vii). 





Furthermore, there is a sufficient margin to allow for any conceivable 
inaccuracy of drawing. 

This critical case clearly renders the first solution inadmissible. No 
number of such particular cases can prove the second solution correct, 
ind that no more complicated transition occurs for any values of the 
parameters. It is conceivable, for example, that instead of a limit cycle 
| lorming through the col on CP (Fig. 2c) and shrinking up round the focus 

as CL is approached with x decreasing, a second-order limit cycle might 
form round the focus and split as x decreases into a stable cycle which 
shrinks up round the focus, and an unstable cycle which expands until it 
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disappears by passing through the col. Any such possibility would be 
extremely difficult either to prove or disprove by numerical methods since 
it would occur when the singular points were very close together, and the 
limit cycles would be very elongated in one direction, the directions of the 
integral curves varying very rapidly near the extremities of the long axis, 
One can only say that a number of other cases have been plotted, and the 
general way in which the integral curve pattern varies with varying para- 
meters appears to be consistent throughout with the second solution. The 
general character of the transitions in the critical region is shown diagram- 
matically in Fig. 5. 
Fig. 6 shows the integral curves with the double singularity for 


F? = 6-625/27 and wz = 00-5166. 


The case F? = 7/27, x = 0-5330 has also been constructed, though not 
reproduced here. It shows the limit cycle forming through the double 
singularity as in Fig. 4, and hence that the value of F, must lie between 
/(7/27) and ,/(9/32). Thus for values of F from zero up to very nearly 
,/(9/32) the final limit cycle arises, as 2 increases from zero, from the 
separatrix when the node and col unite and annul each other. 


Curves of mean-square amplitude 
The general form of the results may be represented by superimposing on 
Fig. 1 the curves of mean-square amplitude for the limit cycles, i.e. 


j= (1/7) { ydt, 


the integral being taken round the limit cycle, and 7’ being the time in 
which the limit cycle is described. These curves are of special interest in 
the application to electrical oscillators since the mean-square amplitude is 
one of the quantities which can easily be measured. 

Such curves have been plotted by Cartwright, and Fig. 4 of her paper (5) 
requires to be modified in the light of the present result. The hypothetical 
curves are shown in Fig. 7, in which only the region near the vertex of the 
ellipse & is shown, the relation to the complete figure being indicated on 
Fig. 1. 

The remainder of the figure is unaltered. 

It is clear that the region between F? = (4/27)+8 and F? = 1/4, in which 
Cartwright’s solution indicated hysteresis phenomena, is eliminated. Over 
the whole range from F = 0 to the value F = F, the synchronized oscilla- 
tion gives way to the combination oscillation represented by the limit cycle, 
as x increases for a fixed F, at the point where the resonance curve meets 
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the upper boundary of the ellipse &, the transition taking place at the same 
point in the reverse direction when x decreases. 
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Fic. 7. Hypothetical curves of mean-square amplitude of limit cycles. 
The resonance curves and curves of mean-square amplitude of the limit 
cycles are shown on a portion of the (x, y)-plane. The resonance curves and 
the ellipse are shown by continuous lines, the curves for the limit cycles 
and the transition locus (T.L.) on which the limit cycle passes through 
the col, are shown by broken lines. The hysteresis effects for increasing and 
decreasing frequency are indicated by lines parallel to the y-axis representing 
the jumps in mean-square amplitude. 


There still remains a small range from F = Fp to F = ,/(8/27) in which 
a hysteresis effect occurs. For decreasing values of x the limit cycle gives 
way to the lower stable oscillation on y = 4, and on the lower boundary 
of € the amplitude suddenly jumps to the upper stable oscillation. For 
increasing values of « the upper stable oscillation persists until the upper 
boundary of & is reached when the amplitude drops abruptly to the limit 
eycle if F? < 9/32, or to the lower stable oscillation if F? > 9/32, the latter 
giving way to the limit cycleon y = 4. These effects are indicated on Fig. 7. 
Conclusions 

In conclusion, I wish to acknowledge my indebtedness to Dr. Cartwright’s 
work, without which the present paper would not have been written. The 
fact that she overlooked the second of the two possibilities which have been 
discussed merely emphasizes the inherent difficulty of dealing with limit 
cycles arising through a col or higher-order singularity. This difficulty has 
been increased by the absence of a suitable pictorial representation which 
would render more apparent the possible ways in which a transition could 
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be achieved, and particularly the requirements of continuity. The repre- 
sentation on the (2, F)-plane with the parameters as coordinates, different 
types of integral curve pattern corresponding to regions of the plane, 
bounded by transition loci, is thought to be a useful contribution in this 
connexion. 

Legarding the application to electrical oscillators, it has been shown that 
hysteresis effects with increasing and decreasing frequency of applied e.m{f. 
are confined to a much narrower frequency interval and are more limited 
in character than formerly appeared to be the case. 
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APPENDIX 
Proof that no limit cycle can occur for F? < 4/27 when three 
singular points exist 

Although this is implied by the work of Cartwright (3), (5), she does not give 
a detailed proof, and no such proof seems to have been published. 

With the equations (1) and (2) in the form assumed here, symmetry about the 
b-axis shows that no limit cycle can occur with x = 0 for any F. In passing, this 
observation, combined with considerations of continuity, suffices to show that for 
any fixed F there exists a finite interval of values of 2 about « = 0 in which no limit 
cycle can occur. For x increasing from zero, a limit cyele might conceivably form 
through the col surrounding all three singular points. This must be described clock- 
wise. 

On the other hand, if a singular point is taken as origin of polar coordinates (R, ®), 
(b),¢9) being the coordinates of the singular point referred to the original axes, we 
obtain 


® = —xr+b?sin 20+ Rb, sind, 
where 6 O—d,. With 0 i7 this gives 
® = —ax+b2+(1/V2)Rby, 


and this is certainly positive if b? > x, i.e. yy > 2, and R > 0. Thus if a radius 1s 
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jrawn from the stable singular point outwards at the angle @ = 47 the integral curves 
ross this radius anti-clockwise, and therefore no limit cycle can exist which surrounds 
ill three singularities, since we always have y) > x for the stable singular point, when 
three singularities exist. This shows also that there cannot be a limit cycle which 
surrounds only the stable singular point. 

To show that there is no limit cycle which can only surround the unstable singu- 
ity, denote the singular points by C, S, and U, and let A be a point on the segment 
fthe positive b,-axis between the two branches of the 6 = 0 curve, i.e. in the region in 

which 6 is positive with OA 1/v3 which is always possible. Draw the circular are 
centre O and radius OA from A to meet the circle d 0 at B, then it is easily seen 
that integral curves cross the boundary of the sector OA B containing U outwards. 
Therefore, any limit cycle surrounding U only must either lie wholly within this sector 
r wholly outside and surrounding it. The first cannot occur since the limit cycle 
would lie wholly within the circle centre O and radius 1/v2. If b, P (b,, 65), 
h Q (b,, b,) Bendixon’s first theorem requires that 


{ [(E+22) aay, 


cb, | éby 
should vanish when taken over the region interior to a limit cycle, and this integrand 
is positive throughout this circle. The second possibility that a limit cycle surrounds 
the sector OA B is impossible because lines from the boundary of the sector to C, 
along the radius OC, and along the arc of the circle 4 0 from B to C, are crossed 
by the integral curves in opposite directions, i.e. in both cases outwards from the 
rea OBC 


Therefore t 


1ere can be no limit cycle. 











ON DIFFUSION FROM A CONTINUOUS POINT SOURCE 
AT GROUND LEVEL INTO A TURBULENT 
ATMOSPHERE 


By D. R. DAVIES (The University, Sheffield) 
[Received 22 August 1952] 


SUMMARY 


A theory is developed to predict the turbulent dispersion into the atmosphere 

a cloud of air-borne particles of matter, emitted uniformly and continuously from 
a fixed point source at ground level. The eddy diffusivity method is employed, and 
the work is an extension to three dimensions of the two-dimensional theory by 
K. L. Calder. The coefficient of lateral diffusivity at a given point is assumed to 
depend on both the vertical and the horizontal components of the distance of the 
point from the central line of the cloud emitted by the source, and an empirical 
method of specifying the functional dependence is suggested. In conditions of 
neutral atmospheric stability the formula obtained for the distribution of concen- 
tration in the important case of aerodynamically rough flow near the surface leads 
to theoretical results concerning cloud height and width which are in very good 
agreement with the observational data. 


1. Introduction 
It is a fundamental problem in the study of turbulence in the atmosphere 
to devise a theory which will enable us to predict the diffusion of a cloud 
of air-borne particles of matter from a fixed, continuous, uniform point 
source. Even in the particular case of homogeneous, isotropic turbulence 
the present state of our knowledge is not sufficiently advanced to allow an 
exact solution to be found. The problem of a true understanding of the 
mechanics of turbulence in a boundary layer is probably the most difficult 
of all and an exact solution is still far removed from our grasp; it is 
this problem with which we are primarily concerned in the lower layers 
of the earth’s atmosphere. However, practical problems can be quite 
successfully approached using the method of eddy diffusivities and assuming 
that transfer of mass is directly related to transfer of momentum in a 
boundary layer. The technique has been fully described by O. G. Sutton (1). 
The application of this method leads to successful prediction of mean 
diffusion properties in terms of the degree of turbulence of the atmosphere 
and the roughness characteristics of the underlying surface (both expressed 
in terms of the mean velocity profile), at least when the atmosphere is in 
a condition of neutral stability. 

By applying G. I. Taylor’s earlier statistical theories of turbulence to the 
case of atmospheric flow, O. G. Sutton (1) has obtained theoretical formulae 


(Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 2 (1954)] 
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for the mean spread of smoke clouds from point and line sources. Sutton 
finds it convenient to define the boundary of the cloud at a given distance 
down-wind as the locus of points at which the concentration falls to one- 
tenth of the peak concentration at that distance. In neutral stability 
conditions over a long-grass surface the agreement with experiment of 
theoretical width and height of cloud given by Sutton at a distance of, say, 
100m. from the source is reasonably good for practical purposes; viz. the 
predicted dimensions of cloud width and height are 40m. and 13m. com- 
pared with the observed values of 35m. and 10m. respectively, these being 
time means over periods of at least 3 minutes. 

A later alternative approach by K. L. Calder (2) to two-dimensional 
problems, such as diffusion from a very long line source placed normally 
to the mean wind direction, has been particularly successful in the case 
of aerodynamically rough flow; the agreement obtained between theory 
and experiment is excellent. In this work the only assumptions involved 
are that (a) mean mass and momentum vertical transfers are described by 
the same diffusion coefficient, and (6) the horizontal Reynolds eddy shearing 
stress is independent of height; no use is made of mixing length or correlation 
concepts. O. G. Sutton (3) has also refined his two-dimensional solution 
with the aid of the concept of macro-viscosity to obtain excellent agreement 
with experimental results. 

The aim of this paper is to consider an extension to three dimensions of 
Calder’s two-dimensional theory, and the first step in the process is to 
determine the functional dependence of lateral diffusivity on the relevant 
variables. The result is then substituted into the diffusion equation and 
a particular solution obtained satisfying the conditions demanded by a 
uniform, continuous fixed point source, situated at ground level on a 
horizontal surface of uniform roughness. 


2. The equation of diffusion incorporating the lateral effect 

The usual approach following the eddy diffusivity method is made. Let 
U denote the mean wind speed at a vertical height z above ground level in 
the positive sense of the x-axis of the rectangular axes Ox, Oy, Oz and 
suppose that K,, K, are the usual coefficients of eddy diffusivity for mass 
transfer in the positive y and z directions. Ignoring the K,, coefficient,} the 
equation satisfied by the mean concentration x(x, y, z), i.e. the weight of air- 
borne particles of matter per unit volume, is 


7OY Cc — ra 7 
l : : = K,- tT K,- ° (2.1) 
Gx =6ey cy) © ez oz 
+ We therefore assume a mean down-wind drift cf the particles in the x direction, with 
eddy diffusion in the y and z directions. 
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Following Calder’s formulation of the two-dimensional problem we write 3, 1] 
© (=)" (2.2) , 
U; 4 
and A, = «, Ui“. (2.3) 


: — : ! 
In these expressions U, is the value of the mean wind speed at some standard | - 
reference height z,; the quantities m, n, and a, are constants for a given 
measurements, and the relation between m and n is m = n/(2—n). The resi 


interpretation of these factors in terms of two-dimensional theory has been of ° 


whi 
turbulent state which may be evaluated numerically from meteorological | anc 
given in detail by K. L. Calder (2). The formulation of the two-dimensional 
problem by O. G. Sutton (3) differs from this and leads to mathematical ~ 
complications in considering an extension to three dimensions; we shall 
therefore not follow Sutton’s formulation in this paper. 
As equations (2.2) and (2.3) represent a mathematical model of diffusion 
which agrees very well with experimental results for two-dimensional 
problems, it appears that a three-dimensional model will best be defined by 


accepting equations (2.2) and (2.3) and formulating an expression for K,. Th 
It has been assumed by the author in a previous paper (4) that K,, may be ph 
represented as a simple power law of z, the height above the surface; in fact 
it was discovered that the law K, x z!~™ satisfied the conditions of the ” 
problem but led to very considerable mathematical difficulties. 

These mathematical difficulties are now easily overcome, as we shall see 
in section 3, by writing 

K, = a, U}-"|y|*z7, (2.4) { ar 

i.e. we now assume that K,, depends on y, the lateral (horizontal) component F 
of distance from the central axis of the cloud emitted by the source and also | 
on 2, the vertical height above the surface: but the dependence of K, on z 
is still assumed to be according to equation (2.3). _ = 

Some tentative physical justification for the form assumed for K,, in 
equation (2.4) will be given in section 4, but at this stage we regard (2.4)as_ | of 


an assumption, governed by mathematical convenience, leading toasuitable | G 


solution of the basic differential equation (2.1). a 
Substitution of equations (2.2), (2.3), and (2.4) into (2.1) leads to the 0! 
equation ; i : y 
zm CX — bl zl-m©X 2... y X27 cx ’ (2.5) 
Ox oz oz} = oy cy ; 
thar = am/JTn aed am/ [Tn 26 
where A =a,27/U}, B= a,zq/U3, (2.6) 


and positive values only of y are to be considered; A is known from previous 
work and B, «, and y are to be determined. 






























ON DIFFUSION FROM A CONTINUOUS POINT SOURCE 171 


rite 3, Diffusion from a continuous fixed point source at ground level 
Following the method introduced in (4), we write 

99 

2.2 y 2B 21 

- = a U, = = v, (3.1) 

lard and x = Car gf(u, v), (3.2) 


ven where the constants £, q, and nm, are to be determined in addition to a, y, 
‘ical | and B. When equations (3.1) and (3.2) are substituted into (2.5), the 


The resulting equation is seen to be separable only if y = m. With this choice 


een of y the equation for ¢(u,v) then becomes 
nal ,  gto-t Cus rani 
i Ne x" 1), u 1 2) 
1Cal 0 ( F ’ 
| cl ct 
hall at 
' A Q2(ayalia)(28 2m—1 3 4(No 2/q) © op | B(B m)(vx! a\B 2m S)/Bopling—rig) OY mT 
Cv? ov 
s10n ¢ ¢ 
C Cw 
mal + B—| (ua! a) aping—2ig) . {3.3) 
1 | cu ou 
5 | 
K This equation is found to be quite tractable; we also note here that some 
“a physical basis for writing y = m is put forward later in section 4. 
7” The indices of x on the right-hand side of equation (3.3) are now arranged 
the to be equal to (n»—1), and thus 2 is cancelled out. We have 
(28—2m—1 2 
Nn,— 1 a ene 4 Ny — (3.4) 
- Bq 4 
» 
} x — ~ 
2.4 and Ny—l = -+n,—-. (3.5) 
q ¢ 
ent . ‘ 5 = . 
: From (3.4) and (3.5) we find that 
LiISO 
q = (2—a) (3.6) 
mn 2 
and B = (14+2m)/(2—a). (3.7) 
jin | A continuous point source is now supposed to be stationed at the origin 
) as of coordinates, emitting Q gm. per sec. of air-borne particles of matter. 


ble | Given that the strength of the source, Q, is constant, the flux of matter 
} across any infinite vertical plane, normal to the x-axis, must be independent 


the of distance down-wind. This continuity consideration leads to the condition 

Oo) | | | w(z) x(x, y,2) dydz = Q for all x values. (3.8) 
r : 0 

»6) | The point source is described by the conditions 

US X¥>0O asx>y>z-), (3.9) 


; and x>0 as 2-0 withy >0,z> 0. (3.10) 
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As in previous work (4), substitution of equations (3.1) and (3.2) into 
(3.8) leads to the condition 


oO x 


dapeettevonsnina | yb(u, v)v™+1-PB dudv = Q. (3.11) 


0 


CU, 


The index of a has to be arranged to be zero, whence 
Ny = —1/q—(1+m)/B¢q. (3.12) 
Substituting the expressions for 8 and q from (3.6) and (3.7) gives 
1/q+n) = —(1+m)/(1+2m), (3.13) 


which results from continuity considerations in the down-wind direction. 
We return now to the equation for x, which becomes, when x has been 


eliminated, 
l Cus Ob 
Ng b—-—(U— +0 — 
qd cu ov 


= A| B2y(28-2m- wpe +88 m)yB-2m— wp 4 Bo yx (3.14) 
cv cu ou 
Writing % = F(u) Gv), (3.15) 
equation (3.14) separates into the two equations 


d2 G dG v dG 


Ap? (2B—2m-1)/B = )B— 2m—1)/B —— aa AG = 0, (3.38 

AB Te | AB(B—m) s + _/ A (3.16) 

and B- : u® fF) 1 ine + (A—n)F = 0, (3.17) 
du du q "du 


where f, g, and n, have now been related by equations (3.6), (3.7), and (3.13), 
and A is a separation constant. 

We find first that a particular solution of (3.16) which may be expected 
ultimately to fit conditions (3.9) and (3.10) is 


G — ede — eb! mae? (3.18) 
if A = —(l1+m)/(1+2m), (3.19) 
and b = 1/AB?q?. (3.20) 


We note from the above analysis that 

Bq = 1+2m, 
and in addition we note that Calder (2), in his two-dimensional work for 
the case of aerodynamically rough flow, writes 


kK, — = MU, ay Mzi—™, (3.21) 
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where M = 22”"/mq’?. The parameter z, defines the surface roughness and 
depends on the mean velocity profile over the anticipated approximate 
height of the diffusion zone (2). 


Comparison of equation (3.21) with (2.3) and (2.6) shows that A = M, 
so that 
b = 1/M(1+2m)?, (3.22) 


s in Calder’s two-dimensional theory, as was to be expected. We are left 
to deal with equation (3.17) for F which, with the aid of substitution 
3.19) and equation (3.13), becomes 


l lF 1F F 
Pd ed PO (3.23) 
du du qdu q 
An exponential solution again exists of the form 
F = e-wBe, (3.24) 


[t is of interest to note here that this exponential solution only exists 
provided the continuity equation (3.13) is valid. A particular solution of 
the diffusion equation (2.5) is thus 


l \ /zA+2m)\ 9 ™ l y? er 
(ara amp)(S)|e*°| -(aal(=)} 29 


where 1/g+n) = —(1+m)/(1+2m). We shall find later that for values of 


xX Cx™ exp 








met with in the atmosphere é/eéy = 0 at y = 0. The constant C is easily 
determined from the continuity equation (3.13) in terms of the parameters 
Band q. The result obtained is 
Qzi" (1+ 2m)[M(1+-2m)?]-A+ma+em 


“= 90, BD : 3.26 
2U, BVP (q>)P[(1+m)/(1+ 2m) jq?44 (3.26) 





The solution (3.25) contains at this stage two undetermined parameters 
Band g. We note also that the net lateral flux of matter across a vertical 
plane parallel to the mean wind direction and distant y from the z-axis is 
given by 


. 


[ [- K, =] dadz. 


. « 


0 z=VU 


This expression, when evaluated by using (2.4), (3.25), and (3.26), is found 
to be independent of y, g, and B, and is equal to }Q as anticipated, the 
result serving as a check on the preceding analysis. Moreover, if we integrate 
equation (3.25) with respect to y from zero to infinity and double the result, 
We arrive at a solution, independent of B and q, for the infinite cross- 
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wind line source problem which is in exact agreement with that found by 


Calder (2). 


4. Evaluation of the turbulence parameters B and q 

Our present understanding of turbulent diffusion processes does not 
allow an exact determination of the factors « and a, in the expression used 
for the coefficient of lateral eddy diffusivity, K,,, and we now assume that 


x = 1—2m,t (4.1) 


1 yp 2] =n) 
and that oo | ; (4.2) 
a. 2 


w* 
where (v’2)? and (w’2)! are the mean eddy velocities over a given convenient 
period of time at a given point in the lateral and vertical directions. The 
mean value of v’2/w’? over the approximate anticipated height of the 
diffusion zone is then calculated, and indicated in equation (4.2) by square 
brackets []. The lateral eddy diffusivity law is consequently written in 
the form K, = a, Uj-*2*y'-*; (4.3) 
the a, coefficient is determined from equation (4.2), the factor a, having 
been computed from Calder’s theoretical work (2). 

Some indirect physical arguments may now be put forward as partial 
justification of the assumed form for K,, given in equations (2.4) and (4.3), 
and in support of the values adopted for «, y, and a,. 

Firstly, it is clear that an important factor to be considered in a con- 
tinuous diffusion process is the distance from zero level of diffusion; thus 
the power law K, = a, U}-"z!-™ has been found to be appropriate in two- 
dimensional problems involving continuous vertical diffusion from an 
evaporating surface. We now take this K, law (2.3) to be a reasonable 
working hypothesis to describe diffusion in the vertical direction from a 
point source. In this case of diffusion from a continuous point source we 
have assumed in equation (4.3) that the diffusivity K,, describing the degree 
of diffusion in the lateral direction, at a given height z, is an increasing 
function of the horizontal component of distance y from the central line 
of the cloud. The physical interpretation one has in mind of this increase 
of K, with y is that progressively larger eddies dominate the outward 
diffusion of each successive group of particles as the time of diffusion of 
each group increases, y being the effective variable in the present context, 
as we are considering the average dispersion over a long period from a 
continuous source. We should also expect the degree of lateral diffusion 
(expressed in terms of the K,, law) to be modified by the z value, in view of 
the physical idea that as the surface is approached so the diffusive action 


+ The parameter m is less than } in neutral stability conditions. 
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for a given y value decreases with increased breaking down of the eddies 
near the surface. 

Secondly, we anticipate intuitively that the mean dimensions and veloci- 
ties of the eddies become increasingly smaller in the vertical direction as 
the surface is approached, the constraining influence of the surface affecting 
the lateral mean dimensions and velocities to a much lesser degree. We thus 
picture the eddies becoming progressively flattened in vertical extent as we 
approach the surface. The available observational evidence (1) in support 
of this notion indicates that, although above 30m. height the turbulent 
energy tends to become equipartitioned in the x, y, z directions (as long as 
we consider only a time scale of less than a few minutes), at low levels it is, 
in most atmospheric conditions, very definitely not equipartitioned, the 
magnitude of the energy of the lateral eddy velocities being markedly greater 
than that of the vertical components. This effect is accounted for in a very 
approximate and tentative manner by equation (4.2) for the ratio a,/a,. 

A partial justification of the values adopted for «, y, and a, may now be 
, and K, given by equations 
2.3) and (2.4) are, of course, dimensionally equivalent, and if we make the 


given. The eddy diffusivity coefficients K 


not unreasonable assumption that the quantities a, and a, also have the same 
dimensions, then it follows immediately from (2.3) and (2.4) that «a = 1—2m 
as assumed in equation (4.1)), since for mathematical convenience, 
we have written y = m. We then have the relation g = 2—a = 1+2m. 
It is also interesting to note that a part of both diffusivity expressions sug- 
gested then increases with height according to the power law z™, which also 
describes the increase of the mean wind velocity with height. Now it is 
known (1) that mean eddy velocity is roughly proportional to mean wind 
velocity and hence the power law represents approximately the increase with 
height of mean eddy velocity. In this very approximate manner the 
process involving the breaking down of the eddies as we approach the 
surface is thus allowed for in the eddy diffusivity expressions. 

On physical grounds it is reasonable to suppose that the ratio of a, to a, 
is roughly proportional to the ratio of the mean eddy velocities (v’?), (w’?)! 
in the corresponding directions. It has in fact been shown by Sutton (5) 
that when the mean wind speed expression is assumed to be independent 
of height on the left-hand side of the diffusion equation, a, is proportional 
to (w’2)!-». We now assume that this result is true in general, and that a, 
is proportional to (v’*)!~" with the same constant of proportionality, and 
equation (4.2) follows. Indirect support for the use of equation (4.2) also 
comes from the fact that the ensuing theoretical results are found to be in 
good agreement with the results of experiment in neutral stability conditions 
over a long-grass surface. 











176 D. R. DAVIES 


5. Comparison of theoretical and experimental results for diffusion 
from a continuous uniform point source at ground level 

A standard set of experimental values} specifying diffusion from a con- 
tinuous fixed point source has been given by Sutton (1). These values are 
quoted in Table 1 together with the theoretical results obtained from the 
analysis given in this paper. Cloud width and height (as already ‘defined 
in section 1) at a specified distance down-wind from the source, say 100 m., 
as suggested by Sutton (1), and peak concentration (i.e. the value of X 
along the x-axis) are taken as the reference standards. Theoretical formulae 
for cloud width and cloud height are easily obtained from equation (3.25): 
if w and h are the width and height of the cloud respectively at a distance z,, 
then w is given by 2[ B(1+ 2m)™@, log, 10]/@+2™ and h is given by 

[.M(1+2m)?x, log, 10}/4+2m), 

Aerodynamically rough flow over a long-grass surface in neutral stability 
conditions is considered; the values of the relevant parameters are given 
in Calder’s paper (2) and shown in the table below. The value of [v’2, w"?}! 
is 1-8, taken as a mean over the lowest 10m. of the atmosphere, as quoted 
by Sutton (3). 


TABLE 1 (a). Numerical values of the physical entities (cm. sec. units) 
employed in the calculation 


Q = 1 g./sec. m = 0-205 %=3 
q’ = 4:42 A = 0:39 U, = 500 
z,= 200 3 \[v2/w?] = 1-8 B = 0-86 


TABLE | (6). Theoretical and experimental results giving cloud width, height, 
and peak concentration from a continuous point source, at ground level, 
in neutral stability conditions over a long-grass surface 


Calculated Observed 
Width of cloud at 100m. _ f ‘ ‘ - ‘ 37-0 m. 35-0 m. 
Height of cloud at 100m... ; : z . ‘ 10-5 m. 10-0 m. 
Peak concentration? from point source, at 100m. down-wind 1-75 mg./m.* 2-0 mg./m.* 
Law of decrease of peak concentration with distance from 
point source. —™ ; ghee St 


An important factor to be considered is the value of n). As we have now 
taken « to be of the form (1— 2m), then the n, value which follows from equa- 


tion (3.13) is given by —ny = (2+m)|(14+-2m). (5.1) 


Using now the value of m (0-205) given by Calder (2) for a distance of 100m. 
travel over a long-grass surface, we arrive at a value for —n, of 1°58; corre- 
sponding to a distance of travel of 1,000m. the relevant m value from 
Calder’s theory is 0-17, leading to a value for —n, of 1-62. The observed 
value, quoted by Sutton (3), for the —n, parameter is 1-76. 

The results obtained are summarized in the table, which indicates that 


+ These refer to average diffusion over periods of at least 3 minutes duration. 
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the agreement between theory and experiment for cloud height and width 
s very good. A slight discrepancy exists between the theoretical and 
served values of the —m, index in the law of decrease of peak concentra- 
tion at a distance of 100m. down-wind from the source; however, the 
wreement obtained between theory and experiment for the magnitude of 
peak concentration at a distance of 100 m. down-wind is good. We note also 
that if we use the observed value 1-76 for the index —n, in order to calculate 
the observed magnitudes of peak concentration at 50m. and 1,000 m. down- 
wind of the source, from the result (quoted by Sutton (1)) at 100m., and 
ipply the appropriate values of m (given by Calder (2)) to calculate theoreti- 
cal peak concentrations at these distances, we find a similar order of agree- 
ment; e.g. at 1,000 m. the theoretical value of peak concentration is 
0-030 mg./m.? compared with the observell value of 0-035 mg./m.* 


6. Application of the power law K,, =- a, U}-"z™y!~2 to the problem 

of evaporation from limited saturated areas 

The important problem of calculating the theoretical distribution of rate 
ofevaporation on an evaporating plane surface, limited in length and width, 
presents great mathematical difficulties. There appear to be two possible 
methods of approach: (a) given a suitable K, 
equation (2.1) satisfying the conditions prescribed on an evaporating sur- 
face, (6); (b) regarding an evaporating area for mathematical convenience 
1s a very close grid of point sources, we seek to evaluate the distribution 


law we need solutions of 


of sources so that the resulting concentration of vapour at all points on the 
evaporating surface is a constant, the saturation value: in method (b) the 
power law K, = a, Uj-"z™y!-*™ can be used for each point source. 
On initial consideration it seems unlikely that the form 

K, = a, U}-*z™y!-* 

will be suitable for the area source problem if we follow method (a), and it 
isnot immediately clear what the form of dependence of K, on y should be 
for these problems. Method (6) leads to a difficult integral equation for the 
source strength distribution function on a limited afea, and this has not 
yet been solved even when the K’s are assumed to be constant. It should 
be noted here, however, that, for reasons described fully in a previous 
paper (4), the law K, = a, U}-"2™ serves as a useful working hypothesis as 
long as we confine our attention to the distribution of rate of evaporation 
on the area and the diffusion in the region near to the surface. The prob- 
lem of evaporation from a limited area bounded by a parabolic curve was 
solved in this way by the author, and the numerical results (based on an 


arbitrary choice of a,) agreed well with experiment in the conditions 
tested (see, e.g., (4)). 
It is also of interest to note here that it is possible to solve, on the basis 


5092.26 


N 
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of the power law K, = a, Uj-"2'-™, cf. (6), the hypothetical problem of 
diffusion of matter from a plane area of infinite length and finite width, 
composed of infinite line sources, the mean wind direction being parallel 
to the sources. The strength distribution of the sources is chosen to be con- 
stant along the mean wind direction and variable across wind, so that the 
concentration of vapour is maintained constant at all points on the area. 
The theoretical dependence of rate of evaporation per unit length of the 
area on the width of area is found to be w!-”, where w is the width, and 
this result has been verified experimentally (6). Too much significance 
‘~annot be claimed, of course, for this agreement, as the theory assumes 
the rate of evaporation to be independent of the down-wind distance from 
the leading edge of the area. The interesting point now is that the analysis 
of the previous paper (6) applies in the case K, = a, U}-"z"y!". The 
point source solution (3.25) is first integrated with respect to x from zero 
to infinity and, using (3.13) and (5.1), we find, for an infinite down-wind 
line source positioned on the x-axis, that y is proportional to 


(l—m)/(1+2m) 
yi roe +2m 
y | 4 sad 


and, along z = 0, x is then seen to be proportional to y”-!. The analysis 
of the previous paper (6) then follows immediately, and we find that the 
theoretical rate of evaporation per unit length of the infinite area is propor- 
tional to w!-™. The analysis can also be made to show that the (1 —m) index 
in the w’-™ power law depends fairly closely on the index of the power of y 
in the K, law. This result thus provides some further support, although 
admittedly of a very indirect and indecisive nature, for the form adopted 
for the K, law. 


7. Conclusions 

The important feature of the investigation described in this paper is that 
a particular solution of the basic diffusion equation (2.1) may be obtained, 
satisfying the conditions demanded by a uniform continuous fixed point 
source at ground level on a horizontal surface of uniform roughness and 
leading to good agreement with the available observational data, if we 


assume that : o_o r 1—n»m,)1—2m 
K, — a, l a, A K, a, l ‘ee aie 


with a,|a, = [v2 /w?}0—™, 
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THE SKIN FRICTION ON INFINITE CYLINDERS 
MOVING PARALLEL TO THEIR LENGTH 


By G. K. BATCHELOR (Trinity College, Cambridge) 


Received 22 August 1952] 


SUMMARY 


The frictional force on a cylinder moving steadily parallel to its length through 
viscous liquid which is initially at rest is determined with reasonable accuracy 
ver the whole range of values of the duration of the motion and for a wide variety 
fshapes of the cylinder cross-section. When the time ¢ is small, the first approxima- 
m gives a force per unit area which is the same as that for a flat plate of infinite 


1e second approximation takes the shape of the cylinder into account and 


the force on unit length of cylinder is determined in terms of the number of corners, 


ind their angles, in the cylinder cross-section ; if there are no corners, the force on 
init length of the cylinder is the same, to this approximation, as that on a circular 
ylinder of the same perimeter. For large values of t the determination of the fric- 


tional force is reducible to that of a potential problem, the solution of which is known 
fora number of different shapes. The approximations for small and large values 


ft for any one cylinder do not overlap but can be joined without much ambiguity. 
For no value of ¢ do the forces on cylinders of different shape (excluding those whose 
wvature is not everywhere inwards) differ by more than about 25 per cent. 


Introduction 
[His paper is concerned with the three-dimensional unidirectional motion 
that is generated by the forced motion of an infinite cylinder parallel to its 
ength in a viscous incompressible fluid. The cylinder and the fluid are at 
rest initially, and the cylinder is suddenly given a velocity W which remains 
steady thereafter; the subsequent motion of the fluid is to be determined. 
\ knowledge of unsteady unidirecticnal flow of this kind may perhaps be 
ised as a guide to the solution of the more difficult non-linear problem of 
the steady flow produced by immersing a semi-infinite body of the same 
ross-section in a uniform stream. The correspondence between the two 
problems was first suggested by Rayleigh, who made use of it for the (two- 
limensional) case of a flat plate of infinite width; the differences between 
nalogous quantities in the two problems for this simple case are now known 
tobe not inappreciable, so that the correspondence must be used cautiously. 
In what follows only the unsteady, accurately unidirectional, flow past 
nfinite cy linders will be considered. 

Several accounts of three-dimensional viscous motions in the above class 
we been published recently. Howarth (1) has considered the case of a flat 
plate of semi-infinite width, the edge being parallel to the direction of 
lotion. Sowerby (2) has found the velocity distribution generated by the 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 2 (1954)] 
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motion, parallel to its edge, of the wedge formed by the intersection of two 
semi-infinite flat plates at an angle (in the sector occupied by the fluid) z/n, 
where n is an integer. Later, Hasimoto (3) and Cooke and Sowerby (4) 
independently obtained the solution to this same problem for an arbitrary 
angle of intersection of the plates, 7+-«,say, Rayleigh’s results corresponding 
to a = 0 and Howarth’s results to « = z. In these flow problems each of 
the two flat sides of the ‘cylinder’ has infinite width, and no length is defined 
by the geometry of the boundary. The region in which the fluid velocity » 
is significantly different from zero continually increases in width (measured 
normal to the boundary), and the velocity distribution is a function of 
x/(vt)? only, where x is the position vector in the plane perpendicular to the 
flow, ¢ is the duration of the motion, and v is the kinematic viscosity. 

A different state of affairs exists when the cylinder is of finite width. 
When ¢ is very small and the boundary layer has a thickness small com- 
pared with a length which is locally characteristic of the cylinder cross- 
section, we may expect that the velocity distribution will be the same as 
if the boundary were a flat plate. However, for slightly larger values of t 
the velocity distribution near any point on the boundary will be influenced 
by the shape of the boundary in the neighbourhood of that point. The flow 
in the neighbourhood of a sharp edge of the boundary has the form described 
by Hasimoto’s and Cooke and Sowerby’s solution, while near parts of the 
boundary where the curvature is finite, we shall see that the flow is the 
same as if the boundary were circular with the same curvature. The com- 
plete solution for a cylinder of circular cross-section has been worked out 
as a problem in heat conduction (Carslaw and Jaeger (5)), and we can use 
it to obtain information about the skin friction on cylinders of arbitrary 
shape, for certain values of t¢. 

For values of ¢ such that the thickness of the region in which w + 0 has 
linear dimensions of the same order as those of the cylinder cross-section, 
the distribution of w is influenced by the shape of the whole cross-section 
of the cylinder, and for these values of ¢ it is clearly difficult to find general 
results. However, for very large values of ¢ the cylinder is immersed in a 
large region in which the velocity has been brought near to its maximum 

value W, and the exact shape of the cylinder clearly has a diminishing 
effect on the velocity distribution as a whole. We shall see that the skin 
friction at these large values of t can be related to that on a circular cylinder 
with the aid of a simple auxiliary potential problem, and again we can make 
use of the known complete solution for a circular cylinder. 

These separate investigations of the cases of ¢ small and ¢ large permit 
a reasonably complete picture of the variation of the skin friction over the 
whole range of values of t to be obtained for a very general class of cylinder 
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SKIN FRICTION ON INFINITE CYLINDERS 181 
shapes. In view of its fundamental place in the whole analysis we begin 


with a description of the complete solution for the case of a circular cylinder. 


A circular cylinder moving parallel to its length 
Let w be the velocity, parallel to the generators of the cylinder, at a point 
in the fluid specified by coordinates 2, y in a plane perpendicular to the 
generators. The equation that w must satisfy, whatever the shape of the 
cylinder cross-section, is 
low cw. Ow 


= om =~» (1) 


v ot Ox*  Oy?* 
ind the boundary conditions are 
a) w = 0 everywhere in the fluid at t = 0; 
(b) w = W, a constant, on the boundary of the cylinder, for all t > 0. 


The equivalent heat-conduction problem is to find the (variable) tem- 
perature distribution in a uniform solid, of thermometric conductivity v, 
which is initially at zero temperature and which is bounded internally by 
1 cylindrical boundary held at constant temperature W. A complete 
solution of this problem for the case of a circular cylinder is described by 
Carslaw and Jaeger (5) (so far as I know, this is the only known complete 
solution for a cylinder whose cross-section involves a finite length),+ and 


is as follows: 


get (Fol kr )¥q(ka) —Jo(ka)¥q(kr)\ dk 


7 JO Seka) Y eka) | 
where 7 is the distance from the centre of the cylinder, a is the cylinder 
radius, and J), Y, are Bessel functions of order zero. The quantity of greatest 
practical interest is the frictional force per unit area (= 7) at the surface 
of the cylinders: 

= CuW/a, say, (3) 


T = p|0w/Or|,—¢, 


where u is the viscosity of the fluid. With the aid of the identity 
J, (z)¥ o(z) —J o(z)¥o(z) = 2/22 
we find from (2) that the dimensionless coefficient C is given by 
Y 4c , k2vt dk 
sf (ud a”) 9 | i seeks 2 Seek (4) 
nm J kl J3(ka)+ Yi(ka) | 
0 

Numerical values of C are shown in Figs. 2 and 3. 

+ Formal solutions for the cases of an elliptic cylinder and a flat strip have been found 


recently by Tranter (6) and Davies (10) repectively, but no numerical values have been 
published. 
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Carslaw and Jaeger point out that the value of C for small values of 
vt/a® (= T) can be obtained from the expansion 
O(T) = (nT)! 4—H(T [a $47 +... (5) 
The first term in this series gives the value found by Rayleigh for the skin 
friction on a flat plate of infinite width, which is what we should expect 
when the thickness of the boundary layer (which is measured by (vt)!) is 
small compared with the cylinder radius. The second term of the series 
(for 7) is the first in which the curvature of the cylinder boundary occurs, 
and gives the second approximation to 7 for small values of T as 
t = pW[(mt)-?+-1/2a]. 
The corresponding second approximation to the total skin friction on unit 
length of the cylinder thus consists of the infinite flat plate value pW (zvt)-! 
integrated over the circumference, together with a constant term zpW. 
We shall see later that the curious result that the addition to the flat plate 
value is independent of the size of the cylinder is not confined to the case 
of a circular section. 
For large values of vt/a? Carslaw and Jaeger give the expansion 
O(T) = 2/[log(47) —2y]—2y/[log(47)—2y]2-+..., (6) 
where y = 0°577... is Euler’s constant. Further information about the 
situation when vt > a? may be obtained by returning to the expression 


(2) for the complete velocity distribution and writing it in the form 





wan Wt 2W je eee dx 


7 , J5(«A)+Y2(KA) Jw’ 
0 
where Kk = k(vt)}, R = r(vt)-, A = a(vt)-!. 
In view of the known expansion 
1~\4 1,)\6 
$oY(2) = (log 42-tyWo(e)+ (He) EM 4 p49) EN 


we can write 


9 
Tole R)¥o(A) —dh(eA)¥(eR) = “log 5+ P(A, xR), 
7 


where the order of P when eA <1 and «R <1 is 

K*A*log(a/r) or «?R*log(a/r), 
whichever is the greater. Hence, since the value of the integral (7) is domi- 
nated by the behaviour of the integrand at values of « of order unity, we have 


4W on r = e-K* dk 
nm a) J2(«A)+Y2(KA) x 
0 


(38) 


j 


w~ W-— 














prov 
the 
be } 
the 


whe 


fi 
that 
com 
not 
(9) | 
the 
tior 
regi 
of t 
was 
eyli 
inc! 


7 


anc 
mo 
cha 
to’ 


wil 


mi 





les of 


Skin 
cpect 
f)?) is 
seTies 


‘Curs, 


(6) 
the 


ssion 


ymi- 


lave 











SKIN FRICTION ON INFINITE CYLINDERS 183 


provided A <1 and & <1. This asymptotic velocity distribution gives 
the correct skin friction at the cylinder at all values of ¢t (which seems to 
be partly accidental; the method of derivation ensured only that it gives 
the correct skin friction when vt > a), so that we can write 


w/W ~ 1—C(T)log(r/a) (9) 


9 9 


when a*<vt and r< vt. 


The meaning of this logarithmic variation of the velocity at points such 
that 7? < vt is evidently that when the region in which w/W is not small 
compared with unity has become large, the velocity distribution at points 
not near the outer edge of this region is approximately steady; the expression 
9) satisfies the equation (1) by making the right side identically zero and 
the left side asymptotically zero (as the expansion (6) shows). Time varia- 
tion of the velocity distribution is more significant in the outer part of the 
region of non-zero velocity, and at points close to the cylinder the magnitude 
of the velocity, but not the dependence on r, changes slowly. This situation 
was to be expected since when the velocity in the neighbourhood of the 
cylinder boundary has been raised to a value close to its maximum, further 
increase in w must occur at a very slow rate. 

The total frictional force on unit length of the cylinder is 

F = 2nar = 2nCyW 
ind in view of (6) we have the result that F > 0 asa—>0O. A rigid line 
moving parallel to itself through a fluid cannot produce a finite rate of 


change of momentum of the fluid, although the total momentum given 
to the fluid by unit length of line after a time ¢, viz. 


t i 
[ F dt = 2xpWa® ( C(T) aT, (10) 
0 0 


T 
will be finite when ¢ is infinite and of order such that a* | (log 7')-! dT’ is 
finite. 


> 


Cylinder of arbitrary shape; vt < I? 

Consider now a cylinder whose cross-sectional shape is arbitrary. As the 
length 7, representative of the cross-section, we shall choose the perimeter 
divided by 2z, since this length is clearly significant for the total frictional 
foree. We shall suppose that the cross-section consists of a number of sharp 
corners separated by lines whose radius of curvature is nowhere of an order 
of magnitude less than /. The thickness (measured normal to the cylinder 
boundary) of the region in which w/W is not small compared with unity is 


measured by (vt)!, and when vt < 1? a boundary-layer type of approximation 
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to the flow equation (1) will be valid. To see this in more detail we transform 
the equation to a suitable coordinate system. 

We define an orthogonal coordinate system (€, 7), in the plane perpendicu- 
lar to the cylinder, such that the curves 7 = const. are all ‘parallel’ (meaning 
that the distance from any one curve to any other along a normal to either 
curve is constant) to the cylinder boundary; a similar procedure has been 
used by Howarth (7) for the more general case of a three-dimensional 
boundary. The element of length in the plane is given by 


ds? = h? dé*+-h2 dy?, 


where h, is constant by definition and can be chosen as unity, and the radius 
of curvature (measured inwards) of the curve » = const. is 


R— (; chy wi 
h, on 


With this coordinate system the equation (1) becomes 
l ew a 1 é (- 7) + 1 @ (i ) 
v ot h, e&\h, c&} © hy en en 
Cw lew 1 é/(1 éw 
needs Nentatnte S a) Peace (11) 
én? Ren cult x) 
and the boundary conditions are 
(a) w= 0 everywhere at t = 0, 
(6) w = W at the boundary (7 = 0, say) for t > 0. 
Now when the region in which w/W is not small compared with unity is 


1 


very thin (i.e. when (vf) 1) we shall have 


c C 


wT 


on C 
and (vt)'(@/én) = O(1), as in standard boundary-layer theory. Hence 
provided h, and éh,/é& are O(1) the first approximation to (11) is 


ip iil (12) 


the solution of which is the distribution appropriate to a flat plate of infinite 


width, viz. w/W = 1—erf dn(vt)-4, (13) 


giving a friction coefficient equal to the first term in the expansion (5). 
The second approximation to (11) under the same conditions is obtained by 
retaining the term which is one order of magnitude (in (vt)!/R) smaller than 


C*w/On?: 


1 Ow Cw 1 Ow 
— em pe 


— BR (14) 
Ry en 
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form { where R,() is the value of R at the boundary 7 = 0. This equation is the 
same as the exact form of the flow equation for the case of a circular cylinder 
licu- | ofradius Ry, and the solution is the relation (2), although only the approxi- 


ning mate form of (2) for vt <a? is relevant here. The corresponding second 
ther :pproximation to the skin friction is obtained by retaining the first two 
been terms in the expansion (5), since the ratio of these two terms is O{(vt)*/ R}, 
on: is , = 
nal and 7 & pW[(avt)4+-4.Ro}]. (15) 


The radius of curvature of the cylinder boundary in general varies around 


the perimeter, so that the second approximation to the total friction on unit 
dius 





length of a cylinder such that h, and éh,/e€ near the boundary are O(1) is 
F Drds = pW 27 (mt)-?+4 > Re* ds|, 


where s represents distance around the boundary and 2z/ is the perimeter. 
Hence the average of the friction around the boundary of cylinders with 
smooth boundaries is 


F/(2al) = (uW/D)[l/(mt)*+-4] (16) 
(11) and is the same, to this order of approximation, as the friction on a circular 


cylinder of radius 1. 

We proceed now to allow for the existence of a number of corners in the 
cylinder cross-section. In the neighbourhood of a corner it is not true that 
h, and ch,/éé are O(1), and (12) and (14) will cease to be valid approximations 
for (11) as soon as (vt)! is of the same order as the (small) radius of curvature 
of the corner. If we suppose that the corner is so sharp that it is possible 
to find values of (vt)! which are small compared with / but large compared 
with the radius of curvature of the corner, then at these values of t the above 
analysis will apply to parts of the boundary not near a corner, while the 
nce corners can be regarded as ideally sharp. The influence of a corner extends 


f order (vt)!, and so for vt < |? the flow near a corner is the 


over a distance u 

, same as if the corner were formed by the intersection of two plane walls. 

12) Hence for details of the flow near any one of the corners of the cylinder we 
can refer to the papers by Hasimoto (3) and Cooke and Sowerby (4). In par- 
ticular we can obtain from these papers the excess of the total skin friction 
13) (per unit length of cylinder), due to any corner of exterior angle 7+-«, over 
5 that for one side of a plane of the same (infinite) wetted area. This excess 
| “s friction, to be written as 


K(a) pW, (17) 


han — ; ; cas ae ; 
where K is a dimensionless coefficient which is a function of « only, repre- 


sents the contribution to the friction on unit length of the cylinder due 
solely to the presence of a discontinuity in slope in the cylinder section. 


14) 
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Hasimoto (3) and Cooke and Sowerby (4) give (consistent) tables of values 
of K(«) for various values of « (in Cooke and Sowerby’s table the values of 
x are rational multiples of 7), and Fig. 1 has been constructed from these 
tables. 

We are now in a position to write down the total skin friction, correct 
to the order (vt/I?)°, on unit length of a cylinder whose cross-section consists 
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Fic. 1. Excess drag due to a corner. 


of several sharp corners (no two of which are within a distance of order (vt)! 
from each other) separated by smooth curves with curvature of order /-1 
or less. The contribution to the total friction from the portion of the 
boundary between two adjacent corners at s = s, and s = 8, is 





82 


2 Se 
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l rds = pW] (52-8) (xvt)'+-3 | RS* as| 


$1 81 


= pW[(s.—s,)/(avt)!+4(0,—4,)], (18) 


a 


where 6,—8, is the angle between the tangents at the two ends of this 
portion of the boundary, while the contribution due to a corner of exterior 
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angle 7-+-a, is n.WK(a,). Hence the total friction on unit length of the 


cylinder is 
F = ¢$7ds pW[ 2a (mt)i+4(27 — Ya,)+ > K(,)], (19) 


— 
n n 


. 


and the average of the friction around the perimeter is 


F/2zxl (uW I)|t (avt)4 4+1/27> (K (o4,)—4a}]. (20) 
n 

Reference to Fig. 1 shows that, when a is positive, K(«)—}a is always 
negative and thence that the friction on cylinders (with inward curvature 
everywhere) with corners is always less than the friction on cylinders with 
smooth boundaries; to turn a corner sharply is more economical on drag 
that to do it gently with the same perimeter.t 

One or two results concerning cylinders whose sections are regular 
polygons, for which, as with any entirely flat-sided figure, } a, = 27, are 

n 
immediately apparent. For a regular polygon with m corners we have 
y = 27/m for each corner, so that 
> K(«,) = m K(2z/m). 
n 

The integral value of m which makes m K(2z7/m) a minimum is 2; its value 
is then 1-0. For m = 3 it is 1-20, and for m = 4 it is 1-24. Thus the flat 
plate of finite width and zero thickness has less drag (when vt < /*) than 
a cylinder with any other regular polygonal section of the same perimeter. 
Indeed, in view of the downward curvature of the curve in Fig. 1 the flat 
plate has less drag than any other entirely flat-sided cylinder with corners 
such that «, > 0.§ To the order (vt/l?)® the average friction per unit area 
on the flat plate of perimeter 2zl is 


F/2nl = (uW/L)[l/(avt)!+ 1/27]. (21) 


+ Which, of course, is far from being the case when the flow is not parallel to the generators 
of the cylinder. 

{ It will be noticed from Fig. 1 that (m/27)K(27/m) approaches a finite limit, equal to 
about 0:25, as m—+>oo. This limit does not give the drag appropriate to a circular cylinder, 
as one might at first expect, for the reason that the formula for the excess drag due to the 
corners ceases to be valid when any two corners are within a distance of order (vt)? from 
each othe e. 

§ The qualification a 0 is necessary because K(«)+ K(—a) is always negative and 
the drag of a flat-sided cylinder could be reduced by adding an angle with outward cur- 
vature ; indeed if the cross-section consists of a tightly rolled line, the drag can be reduced 
to zero. The reason for the large negative values of K(«) as a—> —7 is, of course, that 
when the corner closes up entirely an infinite part of the two planes comprising the corner 
ceases to be wetted area. It is not difficult to establish that as a — —7z, (7+a)K(a) must 
approach a finite limit, which Cooke and Sowerby’s calculations show to be about — 1-37. 
The value of K(«) ceases to apply to a cylinder of finite cross-section when 7+ « becomes 
very small because then the influence of the corner extends over a region with linear 
dimensions « omparable with / 
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Fig. 2 exhibits some of the results that have been established. The curves 
shown all refer to the value of F'/(27uW) under different conditions as a 
function of vt//? and are as follows: (a) C(vt/I?), as given by (4), which is 
an exact expression for the friction coefficient for a circular cylinder at 
all values of vt/I?; (b) 1/(avt)!+-4, which is the value of F (27uW), correct 
to the order (vt//*)®, for any cylinder whose section is without corners: 


40 















































| 
| 
| 
| 
of} 
ss \ 292 approximation 2 7 
\ for cylinder without corners, re ee + 
, 
\ 
7 \ 299 approximation = nt 
for flat Oe 
F me or flat plate, (—b-) = 
2muWw Now arcular cylinder 
x 
fF _ (exact) 
e—— 
oO — | 
| 
ysc approximation 12 = anaes al . 
for all cylinders , (aE c’ /*) 
, | | | 
0 O02 4 Oo 08 oO 12 4 
vtf,2 
Fic. 2. Skin friction on cylinders; vt /l? < 1. 


(c) 1/(avt)'+-1/27, which is the value of F/(2zW) to the order (vt/I?)° for 
a flat plate, and finally (d) //(zvt)!, which is the first approximation to 
F'/(27W) whatever the shape of the cylinder. 


Cylinder of arbitrary shape; vt > /? 

When the thickness of the region of non-zero velocity is comparable 
with /, the velocity distribution near any point on the boundary is not a 
function of the local boundary shape alone and the determination of the 
skin friction is in general a difficult problem. However, the case in which 

t >F is again simple because the cylinder is then immersed in a large 
region of non-zero velocity and the effect of the shape of the cylinder is 
confined to the inner part of this region. The argument is that when 
vt > [?, the velocity at all points in the neighbourhood of the cylinder will 
be close to W, @w/ét will be close to zero, and hence the distribution of w 
in this region will be given approximately by 

sy o2wh 
Cu ew _ 6. (22) 


27 <= 
ox" Oy* 


The inner boundary condition is w = W at the surface of the cylinder. 
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Now the outer boundary condition has circular symmetry so that if the 
region in which (22) is valid extends sufficiently far from the cylinder— 
which is equivalent to requiring that t be large enough—we shall have 
W—w ~ logr (23) 
for large values of r (where r is the distance from some suitable origin 
within the cylinder boundary), provided r is still so small that (22) is valid. 
At values of r at least as small as those for which (23) is valid the distribu- 
tion of w may be visualized as being identical with the distribution of 
height of a soap film which is bounded externally by a ring of large radius 
lying in a plane and internally by a piece of wire, of the same shape as the 
boundary of the cylinder cross-section, which is in another plane parallel 
to the plane of the ring. There is no need to determine in detail the whole 
of this distribution to find the skin friction. The total drag on unit length 


r Ow 
eo —— €8, 
J en 


where ds, dn are elements of the perimeter and the normal to the perimeter 


of the cylinder is 


of the cylinder section. Now since w is a two-dimensional potential function 
it can also be interpreted as the stream-function of a two-dimensional flow 
and it follows that the line integral of the normal derivative of w around 
any curve surrounding the cylinder is independent of the choice of path. 
In particular the path can be chosen as a circle on which (23) is valid. 
Hence the asymptotic law (23) can be written more fully as 


* 
1\—w/W = C(t, l) log(r/b) (24) 
when 7 is large, where the friction on unit length of the cylinder is 


* 
C(t,))uW; 6 is a length which can be called the radius of the ‘equivalent 
circular cylinder’. 

At values of r larger than those for which (22) and (24) are valid the 
velocity varies with ¢ and r in exactly the same manner as for a circular 
cylinder of radius 6, since (24) is now the inner boundary condition for this 
distribution. Consequently the rate of change of momentum of the fluid, 
and hence the total skin friction on the cylinder, is the same as if the 
cylinder were circular with radius 6 and we can write 


on 

C(t,1l) = C(vt/b?). 
To determine the friction on a cylinder of arbitrary shape at time t we need 
only determine the radius of the circular cylinder which gives rise to the 
same velocity distribution in the logarithmic range; the two cylinders have 
the same drag, which is described by (4) or, equivalently, for these large 
values of t, by (6), with the substitution of 6 for a in both expressions. 
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The determination of the constant 6 for an arbitrary cylinder shape is 
equivalent to the determination of the asymptotic form of the stream- 
function ys for two-dimensional circulatory flow around the cylinde: , since 
this function satisfies V7ys = 0, is zero on the cylinder boundary, and tends 
to %& = xlog(r/b) as r/b > «0, where « is the circulation. In this form the 
problem has already been solved for several different shapes of boundary, 
For an ellipse with semi-axes a, and a, we find from the known complete 


solution (8) that as r > « 
2r 
b> « log : 
a,+4a, 


so that b = }(a,+4,). (25) 


The extreme case a, = 0 corresponds to a flat plate of width 2a,; hence the 
drag of an infinite flat plate of width d at (large) time ¢ after it has begun 
to move parallel to itself is the same as that of a circular cylinder of radius 
d/4 at the same instant, and the first approximation to the drag per unit 
length of plate is 2n.W 
= hl ; (26) 


log{64(vt/d?)!—2y 





It is also possible to obtain the value of 6 for a regular polygon, since the 
conformal transformation that will convert this figure into a circle is known 
(and has been used, for instance, by Morris (9), to analyse problems of two- 
dimensional potential flow about a polygonal cylinder). The radius of the 
circle into which an n-sided regular polygon of perimeter 27/ is converted 
by a transformation which leaves the plane at 7 = oo unchanged is found 
to be o _ 

b = nt | | (2sin g)2” dd _ yp mds t/n) (27) 


92inT (li 
. 22/"1'(4+-1/n) 


For n = 2, corresponding to a flat plate, we have b/l = 47 = 0-785, as 


already found above in a different way, and forn = 3,4, we have b// = 0-885, 





0-925, respectively. When n = ©, corresponding to a circle (and the 
limiting procedure is valid here), we have b/! = 1, and since 6// is a mono- 
tonic function of n there is very little variation of the ratio b// for regular 
polygons. The value of b/I for an elliptic cylinder also varies monotonically 
between the values 0-785 for a flat plate and 1-0 for a circle, and it is probable 
that for all cylinder sections not too different from a regular figure the value 
of b// lies between these same limits. 

Fig. 3 shows the variation of C(vt/I?) (i.e. of the friction coefficient for a 
circular cylinder of radius /) as a function of log(vt/l?). The curves repre- 
senting the friction coefficients for cylinders of different shape but with the 
same perimeter 2z7/ are obtained by translating this curve in the direction 
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of the abscissa through a distance equal to log(b?//?) in each case. Thus the 
curve for a flat plate of perimeter 27/ is obtained by translating the curve 
for a circular cylinder through a distance —0-48, as shown in the figure. 
The corresponding curve for any other regular polygonal section lies 
between the curves for a circular cylinder and a flat plate. 

Although no general theory concerning the skin friction at values of vt 
of the order of /* is available, it is clear that the approximations established 
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for vt < /? and vt > [? put narrow limits on the possible values of F'/(27.W) 
at intermediate values of vt; moreover, the exact expression for F/(27pW) for 
1 circular cylinder provides some guidance to the behaviour of F'/(27u.W) 
for any other cylinder. In the case of the flat plate, forward extrapolation 
of the curve in Fig. 2 can be made consistent with backward extrapolation 
of the curve in Fig. 3 provided the value of F'/(2zpW) at /?/vt = 1 is about 
75, It seems that the maximum difference between the values of F'/(27yW) 
for a circular cylinder and for a flat plate—which are likely to represent 
extreme cases, if we exclude cylinders which do not have inward curvature 
everywhere—is about 25 per cent. of C(vt/l?) and occurs when vt//? is of 
order unity. It is curious that the drag of a circular cylinder is at all times 
greater than that of a flat plate of the same perimeter, and in fact greater 


than that of any cylinder whose section is a regular polygon. 
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\OTE ON THE PAPER ‘TWO-DIMENSIONAL PLASTIC 

STRESS SYSTEMS WITH ISOMETRIC PRINCIPAL 

STRESS TRAJECTORIES’+ 
By A. VAN TUYL 
(U.S. Naval Ordnance Laboratory, White Oak, Silver Spring, Md.) 
[Received 28 August 1952] 

A COMPLETE discussion of Case II, pages 7-9, leads to the following result: 
If A(t), F(t), and G(t) are such that H(x,y) = A(u)+ F(x)—G(y) is har- 
monic, where w(2, y) F(x)+G(y), while D(x, y) = (F’+@")/(F?+@") 
isnot a function of w, then A’(t) is a constant. The proof found originally 
was not given in full, a part being omitted because of its length. The 
purpose of the present note is to give the following much shorter proof. 

The proof rests on the following two results: 

(1) The harmonic functions 

H,,, = F'G'A"(u) and H, H, = F'G'[A’(u)—1] 
are both of the form FG’ B(u). 

(2) If F’G’ B,(u) and F’G’ B,(u) are both harmonic, with F and G such 
that D(x, y) is not a function of uw, then B, is a constant multiple of B,. 
This is an immediate result of the relation 

B, B,— B, B, = —3D(zx, y)(B, B,—B, B,) 
obtained from B,V?F’G’ B,— B,V?F’G' B, = 0), since if D(x,y) is not a 
function of «, both sides must vanish. 
It follows from (1) and (2) that A(t) satisfies the differential equation 
A” ld { A’+1) 
— eee | nel es, 
1—A® ~ 3at| "a—If . 


where c is constant. If c = 0, we find A’(t) = constant; if c + 0, 
J ; 
A(t) log sinh e(t—t,)+-constant. 
- 


Substituting the latter, we find that H(x,y) reduces to 
a ea” eee 

log| ec?#—) — e-c24—) + constant. 
- j 
All harmonic functions of this form are known,{ and are such that D(z, y) 
isa function of wu. Hence, A’(t) = constant is the only possibility. 

+ P F. Neményi and A. Van Tuyl, ‘Two-dimensional plastic stress systems with iso- 
metric principal stress trajectories’, Quart. J. Mech. and Applied Math. 5 (1952), 1-11. 

t P. F. Neményi, ‘Recent developments in inverse and semi-inverse methods in the 
mechanics of continua’, Advances in Applied Mechanics II (Academic Press, 1951), p. 144. 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 2 (1954)] 
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A CONTRIBUTION TO THE THEORY OF 
RIGHT-ANGLED JUNCTIONS IN WAVE GUIDES 


By J. D. PEARSON 
(Department of Mathematics, The University, Manchester) 


[Received 9 September 1952] 


SUMMARY 


The propagation of waves through a junction of two rectangular wave guides at 
right angles is considered for the case when the dimensions of one of the guides are 
such as to allow the propagation of more than one mode. Numerical results are given 
when one guide is of square cross-section. 


1. Introduction 
In this paper we consider the mathematical theory of the junction of two 
rectangular wave guides at right angles, the cross-sections of the guides 
having the same width parallel to their line of intersection, as shown in Fig. 1. 
No other restrictions are imposed on the dimensions of the two guides; thus 
the case of more than one propagated mode may be considered. Energy is 
fed in along one arm of guide 1, the other being closed by a reflector. 

Results have been given previously for the right-handed junction of two 
rectangular guides of the same cross-section by Allanson, Cooper, and 
Cowling (1). 

In what follows, the junction of a square and a rectangular guide is con- 
sidered in some detail, since it has properties similar to those of a junction 
of a circular and a rectangular guide, which is a problem of some practical 
significance. In a circular guide, with a suitable choice of frequency, the 
propagated field consists of two modes. The mode with the higher critical 
frequency possesses rotational symmetry, while the lower mode does not 
have this property. In practice it is required that the maximum amount 
of energy shall be transmitted in the symmetric mode, the non-symmetric 
mode being removed by means of a filter. A phenomenon similar to the 
existence of a symmetric and a non-symmetric mode in a circular guide is 
shown to exist in a guide of square cross-section. 

The junction of the two rectangular guides is as shown in Fig. 1. The 
guides are assumed to have perfectly conducting walls, guide 1 being 
terminated by a perfectly conducting plate. Only periodic phenomena 
varying as exp(twt) are considered; for simplicity this factor will, however, 
be omitted throughout. The problem considered consists in determining 
(Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 2 (1954)] 
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the field in both guides when the source of energy is in guide 1 and is situated 
it a large distance from the junction. The source of energy is such that the 
wave incident on the junction in guide 1 contains only the H,, mode, the 
frequency w being chosen to allow the propagation of this mode. 

Two sets of axes O(x, y,z) and O'(x, n, C) are chosen as indicated in Fig. 1. 
The coordinates y, z are related to », f by the equations 


y b—L, 


& = &9—%)- 


The cross-section dimensions of guide 1 are taken to be a by 6 and those of 
guide 2, a by d. 
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It is shown that if H,, is known everywhere the total field can be deduced, 
The method used to determine H,, consists in assuming an expansion in 
terms of principal modes for @H,,/éy in guide 2 and in particular over the 
aperture € = 5b. From this an expression for H, is obtained in the whole of 
guide 1, as a sum of eigenfunctions of a normal cross-section of guide 1. 
Hence expressions for 2H,,/éy at the planes z = z, and z = z, can be deduced, 
and from these an expansion in terms of modes for H, in guide 2 is again 
lound. Thus two equivalent series expansions for the field component H, 
in guide 2 are obtained and an infinite set of equations is deduced by 
equating the coefficient of each mode in these expansions. 

The analysis given in section 3 is valid even if the field in guide 2 contains 
waves in both directions. The particular case when guide 2 is matched is 
considered in section 4. In section 5 the effect of placing a filter in guide 2 
is discussed. In both sections 4 and 5 numerical results are given for the 
ease when guide 2 is of square cross-section. 
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2. Reduction of the field equations 
Maxwell’s equations in M.K.S. units may be written 





oH, OL. 
—?— —_t = —piwH_Ic, 
ey C2 pie xr Cc (1 (i)) 
OE. cE 

2 — — pin : ii 
Cz Cx eae (1 (i) 
ne a (1 (iii) 
ex ey 
SS... = eiwLE,/c, (2 (i)) 
oy Oz 
OH, OH, ' 
geil ae duane ae E “ Z ii 
G2 ox (2 (i) 
ook me clei, (2 (iii) 
ex ey 


where yp is the magnetic permeability and ¢ the dielectric constant. 

The boundary conditions on the walls of the wave guides, for any com- 
ponent of the electromagnetic field, are that either the component itself 
(viz. EZ tangential and H normal) or its normal derivative (viz. E normal 
and H tangential) shall be zero. It follows that the components of E and H 
may be expanded in one of the following forms: 


x 


pa A,,(y, z)eos oo | 


0 ™ 
n= ¢ 
. (3) 
~ NTL 
2 B,,(y, z)sin — 
a 
n=1 


As the discontinuity at the junction of the wave guides is parallel to the 
x-axis, it follows that each term in an expansion of type (3) may be 
considered separately. Equations (1 (ii)) and (2(iii)) then reduce to two 
simultaneous equations for H,, and E,, and equations (1 (iii)) and (2 (ii)) 
reduce to two simultaneous equations for H, and E,, assuming that H,, and 
E,, are known everywhere. 

From equations (1 (i)), (2 (ii)), (2 (iii)) it is seen that the nth coefficient 
in the expansion for 7, is a solution of the equation 


a2 2 2,2 
nu w nr 
24. ( —"F u, = 0 (4) 


Oy? — a2? c* a 


By) 
Oru 





with éu,,/év = 0 on the guide walls; from (2 (i)), (1 (ii)), (1 (iii)) it is seen 
that the nth coefficient in the expansion for E, is a solution of equation (4), 
with uw, = 0 on the guide walls. Thus solutions with H, zero and with L, 
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zero may be considered separately. The incident field is in the Hj, mode, 
which possesses no Z,, component; we therefore take E, zero, and the prob- 
lem is reduced to finding H,,, a solution of the scalar wave equation (4) with 
n= 1, boundary conditions éu/év = 0 on the guide walls, and a given 
incident field in guide 1. 

If d (see Fig. 1) is taken sufficiently large, the propagated field in guide 2 
will consist of three modes, E,), £,,, and H,,, but the combination of £,, 
and H,, will be such as to have no £,, component. This combination of the 
B,, and H,, modes for a guide of square cross-section is the counterpart of 
the symmetric £,, mode of a circular guide. 


3. Equations for H 


The x-component of the magnetic field may be written in the form 


- WH ~ 
H, = u(y, z)sin—, (5) 
. a 
O*u . mu w? 7 
where = at 7-3 u= v. (6) 
Cy? © ez? ct og 


On the walls 7» = 0, » = d of guide 2 
cu/en = 9; 


hence, in guide 2, « may be expanded in the form 


u = > b,(2)x_(n), (7) 
q=0 
where x,(n) = (2/d)*p, cos(qzn/d) 
and Po l (¢q +9), Po = 1/2, (8) 
d 
b, | ux,(n) dy. 
0 


Similarly in guide 1 we may write 


u = > a,(z)h,(y), (9) 
where ub .(y) (2/b)*p, cos = 
y) 
: 
and a,(z) = | uly, z)b,(y) dy. (10) 
0 


From equations (6), (9), and (10) it is easily shown that 


As __)2q, = A, (11) 


C2 
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S272 2 w? 

where AS == ——— 
- ~§2 er: c 

and A, = [—éu/eyy,(y)| 


The solution of equation (11) is 


b 


y=0° 


—,2 


a, = ers [ e-A=(hA, r,) dz—e-Ast 
1 Ye 
where y, and y, are constants to be determined. 

Using the fact that for z > z, all waves must be travelling in the positive 
z-direction, except for the incident field emitted by the source of energy, 
and also the boundary condition éu/éz = 0 at the plane z = 0, it follows 
that equation (12) may be written 


e?(4A./A,) dz, (12) 


21 rag 
- 


a, = ed [ e-N2(hA,/A,) dz—e-' | e-*#(4.4,/A,) de—e-Ne Jewd 4 ./A,) dz+ 
1 Zo 
+By(e'¥ +e), (13) 
the incident field being taken to have an H, component of unit magnitude, 
The coefficients b,(£) of equation (7) may be shown by a method similar 
to that used in obtaining equation (11) to satisfy the equation 


0b, ss 
ae Haba = By (14) 
c n=d 
where B, = |- = xa(n)| (15) 
on wit 
n 
and w= a 


For ¢ > b, b, may be written 


b, = a, eMaS-+- 8, e-Has, 
and it is easily shown that b 
B, 
B,—% = — [ ~~ cosh(11, £) dé. (16) 
a qd 
0 


From equations (9) and (15) we have 


Be) = 21), > do—()_ —(E)_ I. 


a \ éz ez 
8 a 


Substituting this expression in equation (16), it is seen that 


2 sinh(y, 5) (Ca, ea, 2a 
ne or ui(bd)? Pa 2, {(- y() —(3) Je (1 + 
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From (13) it can be shown that 





Ca, a . . ° 
(— 7 *) -— (‘ :) = 28.9 Arol (—)@ sinh Ayo 2)—sinh Ajy 2, |+ 
sis we b,,( ss Pm 1'(8,m,q) 
} 5“ m m m 18 
aay" nA, } 3 24+ ma? /d? » (18) 


where 
T(s,m,q) = [—sinh A, zo(—e7As%1+ (— )me-As20)(— )a— 
—e*1(sinh A, 2; —(—)” sinh A, zy). 


From (17) and (18), and putting 8, e-#= = B), and a,, e#m’ = we have 


finally 


Xm 


* Xg 


p,| ] ¢ 214g) = al l ¢ 2Hyh) 
2 TS (Ps)?Pm mAs Ls, IB %m) 
bd H1 Hy (8?77®/b® + 5) (mx? /d?+-A2) 


mm 





9\1 l 
«\2 : —— «7 ( 
-( ) Aro—=-[(—)* sinh Ayo 2—sinhA,o2z,}. (19) 

ii lq 

Equation (19) represents an infinite set of linear simultaneous equations 
| 1 

containing the two infinite sets of unknown coefficients a,, Bj. The co- 
efficients "a 
square cross-section. ‘wo particular forms of «, will be considered below; 


arise from waves incident on the junction along the guide of 


these are a, = 0, which corresponds to guide 2 being matched, and 


— ’ 9 —2po(p—b) 
xq = 809 Pa é ue? > 


qd 
which corresponds to guide 2 having a filter at the plane £ = p which 
completely reflects the 1,0 mode and leaves the 1, 1 modes unchanged. 


4. Guide 2 matched 
The simplest form of equation (19) occurs when guide 2 is matched, so 
that all waves in it are travelling away from the junction. In this case 


we have 


and equation (19) simplifies to the form 





— bg 2 K* +S (Ps )? Pm( o a Py (m?—+y? *)Dm T (s,m, q) 
(l—e 7) a ar at hot (8? +. R*q?— B?)(s? + m* R* — B?) 
2!(B/r) R-4(q?—y?)—| (— sin azy—sinaz,|, (20) 
where \ —1Ajo; b,(b) = 6, 


B = ab/z, y=ad/7, and R= Dd/d. 
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The function 7'(s,m,q) may be simplified for particular values z,. If tion 
% = 0, T(s,m,q) reduces to miss 
T(s,m,q) = 3(1—e-*s4), 
and in this case the 7'-junction is reduced to a corner. If z (see Fig. 1) is 
sufficiently large for terms of Ofexp(—A,d)}, s 4 0, to be negligible, then 
T'(s,m,q) reduces to 0, if m is odd and q even, or m even and qg odd, s + 0; 

to —1 if m and q are odd or m and q even, s + 0; to 
—sin(y7)sin(2az)+y7)—? sin(yz)cos(2az)+y7) 


for m even, g odd, or g even, m odd, s = 0; to 





{cos(2az)+y7)— 1}{cos(y7)+ 1}+-i[sin(yz)— {1+ cos(yz)}sin(2az)+-yz)] 
for m and q odd, s + 0; to 


{cos(2az + ym)+ 1}{cos(yz)—1}+-i[ —sin(yz) +{1—cos(yz)}sin(2az)+yz)] 


for m and q even, s + 0. 
It is seen from equation (20) that when the end plate is sufficiently far 


from the corner, so that terms of Ofexp(—A,d)}, s + 0, can be neglected, ” 

the coefficients b, are periodic with period 2z/a with variation in position 

of the end plate. It is also seen that |b,| is periodic with period z/«; and} S 

since the energy in any one mode is proportional to |b?|, it follows that “ 
as 


the transmission coefficient of a propagated mode is periodic with period 
equal to the wave-length of the incident field. ™ 
Approximate solutions of equation (20) were computed for the case when 
terms of Ofexp(—A,d)}, s ~ 0, could be neglected, for various positions of 
the end plate over one-half wave-length of the incident field. The summa- 
tion over m was approximated to by retaining only the first five terms, so 
that the infinite set of equations was reduced to ten equations. The sum- 
mation over s was approximated to by the first fifteen terms, and Newton’s B 
integral approximation was used to estimate the remainder. p 
The values of the parameters in equation (20) used to compute values 
of b, were as follows: 


} 

R = 0-3333, 6 =0-46413,  y = 1-93878, B=1. | 
These correspond to guide 1 having dimensions a = 3 in., b = 1 in., and g 
guide 2 dimensions a = 3 in., d = 3 in., with an incident field of unit 1 
amplitude of wave-length 3-5 in. (8-89 centimetres). 


Values of |b,| are given in Table 1 and values of the energy carried by 
the propagated modes in guide 2 are given in Table 2. 

To check the accuracy of the computation, the coefficient a, of the | 
equation (9) was determined using the calculated values of b,,, and from a) | 
the reflection coefficient was obtained. The agreement between the reflec- 





1) is 


then 


ues 


nd 


nit 





RIGHT-ANGLED JUNCTIONS IN WAVE GUIDES 201 


tion coefficient calculated in this way, and that calculated from the trans- 
mission coefficient was in all cases within 2 per cent. 


TABLE 1 























| in | , | 
0 ) b,| | by Og 
IO 547 “402 617 037 
I*2 457 17048 333 } "105 
79 I*104 257 "103 
I 744 *706 544 066 
mn ¢ 
TABLE 2 
0/4 Po Py 
1'O °717 *II3 
I°2 *209 763 
I*4 "143 849 
6 | *553 *347 





5. Filter in guide 2 

In this section we discuss the effect of placing a filter in guide 2 at a 
sufficient distance from the junction for the evanescent modes arising 
at the junction and the filter not to interact appreciably. The filter is 
assumed to be such that the £,, mode is completely reflected and the 1, 1 
modes transmitted unchanged. 


If the filter is placed in guide 2 at the plane £ = p, we must have 


9 
, »—2t9 
XQ By e- HP 


or Yo - Boé 2Hol(p b). 


and substituting in equation (19) we obtain 


B e—*Holp—0 
‘ 0 _ ole 
) } qo 9 1 
p, \(1—e?#e") * (1—e-2H0?) 
1 


jon oy (p. )*p (s* B?)3(m y*)! , 
eo (s,m, +8, Boe 
R? “q" 87) / (s?-+ m2 2 R2_ —p? 4 m Q[Bm Om Boé 


—2po(p—b) | 


A L, ($4 

2! (8/7) R-4(q?—y?)-[ (—)@sin azy—sinaz,|. (21) 
Solutions of equation (21) were computed by methods similar to those used 
to obtain solutions of (20), for various positions of the filter in guide 2, 
and the end plate in guide 1. It is seen that the energy transmitted in the 
1,1 modes is periodic for variation in position of the filter, with period z/« 
The amount of energy in the 1, 1 modes, for different positions of the filter, 
and a fixed position of the end plate, was found to vary about the value 


of the energy in the 1,1 mode for the matched guide with the same end- 
plate position. 
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Values of the energy transmitted in the 1,1 modes for various positions 
of the filter and end plate are given in Table 3; the values of the parameters 
taken in computation were the same as those used in calculating Tables | 
and 2. 




















TABLE 3 
l 
¢ 
Zo/d Bt ° | Air 7 an 
1'O *223 | “128 "007 "099 
I°2 672 735 | "923 762 
I*4 *742 | “S31 "965 | °857 
1°6 °27 242 5 *297 
74 4 570° a) for 
[id = —2uo(p—b).] | de 
P,, transmitted. 
; su 
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SUMMARY 

A method is developed for determining the pressure distribution to second order 
for a certain class of wings which require the use of three independent variables to 
describe the flow. 

The flow investigated is that of an inviscid perfect gas which is everywhere 
supersonic. 

The method is used to obtain the pressure distribution to second order for a 
trapezoidal wing which has the ends raked outside the Mach cones from the tips 
and which is rolling at constant angular velocity. The rolling moment coefficient 
is also found. The method also yields the position of the shock surface, to the first 
order in the disturbance velocity, in the case of compressive flow and the boundary 
of the region of expansion in the case of expansive flow. 

The mathematical complexities encountered in applying the method to more 
complicated wings are exhibited by a discussion of a rectangular wing. 

1. Introduction 

Up to the present the main non-linear studies of three-dimensional wing 
theory in supersonic flow have been those made on conical flows. The prin- 
cipal results were obtained by Lighthill (1), Moore (2), and Tan (3). 

In the work that follows a method is developed which utilizes Lighthill’s 
technique for ‘rendering approximate solutions to physical problems uni- 
formly valid’ (4) and which is applicable to a certain class of wings that 
require the use of three independent variables to describe the flow. 

The method is used to obtain the pressure distribution to the second order 
for one case, and it will be shown that the practical limitations of the method 
depend on whether or not the first-order solution can be expressed in rather 
simple mathematical terms. 

2. Equations of motion for steady irrotational flow 

Cartesian coordinates Z, 7, Z are introduced with the origin at the foremost 
point of the wing and with the axes orientated so that the uniform undis- 
turbed flow, of velocity UV and Mach number M, is parallel to the Z-axis, 
and so that the Z-axis is perpendicular to the plane of the wing. It is assumed 
that the flow is irrotational, which implies the existence of a velocity 
potential Q. It is convenient to work with a normalized perturbation 

otential defined by 1217 
Q = U(é+9). (1) 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 2 (1954)] 
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The equation governing the motion of the fluid may be written in terms 


of & as (5) 


where B = (M?—1)! 


and y is the ratio of specific heats at constant pressure and at constant 
volume. 


3. Boundary conditions 
Consider first a wing whose leading edge is the j-axis and whose upper 
and lower surfaces are given for all and 7 by 


Z = ¢, S,(@,9)H(z) (3) 
and z= ¢,S(%,9)H(z), (4) 
where H(z) is the Heaviside unit function defined by 


and where, for all 7, 
S,(%,9) and S(%,7) are analytic functions of and 7 for # > 0; 


e8,(0, 9) 


Cx Cx 


a > 


i me 
' = pny as : eS,(0, 7) 
S,(0,97) = S(0,9) = 9; —— 
e, S, > «8, and e, and e, are real numbers. The subscripts u and / refer 
to the upper and lower surfaces respectively. 

We seek a solution to equation (2) which satisfies the following conditions: 


(i) No propagation of disturbance upstream. This condition implies 
that (%, 9, Z), = o(@, 9,2), = 0, F < 0. 

(ii) The flow is to be tangent to the surface of the wing. Under the 
assumption of analyticity of S,, and S,, this condition will be referred 
to the #j-plane by means of a power series expansion in e,, for the 
upper surface and in ¢, for the lower surface. 

From condition (i) it can be deduced that the wing divides the flow into 
two parts, one above and one below the #2-plane, and that the flows in the 
two parts are independent. Thus, two more conditions may be stated as 

(iii) £(%, 7, 2), is to be continuous for all # and 7, Z > 0; J(%, 7, 2), is to 
be continuous for all @ and 7, 7 < 0. 





Se EEnnnnEnEEEEEEEEEee oa 








the 





oTrms 


tant 


he 
ed 


he 


ito 





THREE-DIMENSIONAL WINGS IN SUPERSONIC FLOW 205 


(iv) J(@,9,2), and (%, 7,2), are to be of class C® in regions R,, and R, 
respectively, where R,, consists of the points on and above the Z9- 
plane less the points lying on a surface .@, which passes through the 
g-axis, and R, consists of the points on and below the #-plane 
less the points lying on a surface 4, which also passes through the 
j-axis. (A function is said to be of class C™ in a region R if, in R, 
the derivatives of the nth order are continuous in each of the 
variables.) 

For the present, only the upper half of the flow will be considered and 

the subscripts w and / will not be used. 








4, Scale transformation 
It is convenient to introduce the scale transformation 


i = Bx 
sick iam (6) 


| Under this transformation equation (2) becomes 


My Wz ig, ao\ (eb 
be — Pyy— Ven + Se | —(24- + $4 GAYS +: by +2) + 
Voy /V Wrz B? | 9 ob. ' Bp? ' By v2)( is vy bes) ! 
Qu, us us cr ap? 


B2 — ra gi T Pyy pi a i be. pet apy. py b. - 2.2 v,{l + i) + 


9 Yz\] _ - 

+%bevth(1+33)| a= ©. (7) 

} The boundary condition of tangential flow may be stated as q.n = 0, 

where q is the velocity vector of the fluid at a point P on the surface of the 

wing and n is a vector normal to the surface at P. In terms of the a, y, z 
variables, this condition may be stated as 


[eS (Bx, y), «S,(Ba, y), 1). [B+ by ve 0, (8) 


p* B al. eS(pz, y) 
where S; denotes the partial derivative of S with respect to the first variable 
and[ |.[ | denotes the scalar product of the two vectors. 


5. Reduction to a sequence of linear problems by expansion in 
series 
Following the work of Lighthill (1 and 4), a solution to equation (7) is 
assumed to exist of the form 
fb = eD(u, y,z)+e*h(u, y, z) +e8h,(u, y,z)+... 
x Uutex,(U, Y,z)+e%x,(U, y,z)+... 


, (9) 
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where each series is assumed to have a radius of convergence different from 
zero for u, y, z in R and the 2; are assumed to have derivatives of all 
orders for u > 0, z > 0, and all y. 

Substituting equations (9) into equations (7) and (8), and setting coeffi- 
cients of like powers of « equal to zero, reduces the non-linear problem to 
a sequence of linear problems. The first two problems corresponding to « 
and ¢«? terms will be called the first-order problem and the second-order 
problem respectively. The first-order problem may be stated as 


®,,,,— 9,,— 9... = § ) 
D(u,y,z) = 0, u<o 
@.(u,y,0) = BS(Bu,y), u>OoO ' (10) 


® continuous for all uw, y, and z > 0, and 
® of class C® in R 


and the second-order problem is 


Cx Cx . 
pag: 1 9 1 
Pua dyy— Pez — = e —M ®,,,—2 ey 9,,— 
Dae a2 22 a2 
TW tke ede Mas 
a ™ Ou® dy? az? 





alas = (»— ] + aa ) Duy o.+ 20,,. ®,+ 20, ®,| 


p2 BP 
d(u,y,z)= 0, u<od0 + . (11) 
] , 
b.(u,y,0) = F ©, S, +82, Sy +0, S,—0,, s| , u>0 
2=0 


¢ continuous for all uw, y, and z > 0 
¢ of class C®) in R 

x,(0,y,z) = 0, and 

2,(u, y, 2) to be determined 





6. Solution of the non-homogeneous wave equation with given 
boundary conditions 


If x,(u, y,z) were known, the second-order problem could be written as 


Pua —Pyy—Prus = flu, Y, z) 
$(U,y,2) = g(u,y), u>O 
d(u,y,z)= 0, u<o0 
¢ continuous for all w, y, and z > 0, and 
¢ of class C® in R, 


where f and g are known functions of class C® in R. 
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This problem will be solved by means of the Riesz operator (9 and 10). 
Let V denote a volume which can be decomposed into a finite number 

of simple volumes, V;,i = 1, 2,...,. Let p and q denote functions of wu, y, z 

which are of class C® in each V;, Then Gauss’: divergence theorem may be 


stated 


{ | { [pL(q)—qL(p)] dv = - > IJ (pVq—qVp).v; dS;, (13) 


where, if %; (l;,m;,n;) is the outward unit normal to the surface S,, 
y, = (l,, —m;, —n,) is the outward unit conormal to S;; L denotes the 


: c* Co? o : 
d’Alembertian operator | = ) and V denotes the gradient 
C 


Ee C 7? C cs 
Cc Cc Cc 
operator =» ; ‘ 


0g on oc 


If we let g = G@-D2 = [(u—£)?—(y—n)?—(z—Z)?|@-9?, ~« complex, 


rea > 5, 
Pp $(€, 7, ¢), 


V = volume bounded by the forecone from P(u, y,z) given by 
(u—£€)?—(y—n)?—(z—)? = 0 and the plane é = a,a > u, 


where £, », £ are the so-called ‘running coordinates’, then, since 
L(G&-DI2) = a(a—1)Ge-32 
and L$) =fE0,0), 


equation (13) may be written as 
|} } bE, 0, C)@e-9? dv 
Jd. 
n 
{ | { Go-W2F(E, n, L) dv + > [| (AV Ge-D2— Ge-v2V4).v, dS;. 
t=15; 
Multiplying the above equation through by {27I'(a+1)}-!, where I'(a) is 
the gamma function, yields, in Riesz’s notation, 


I*d(u, y, 2) 
of l ~ rar Y 9 19 Y 
[*+2f(u, y,z)+ — _— S , J (P6VGe-D2_— Ge-D2Vd).v;dS;, (14) 


where I*d(u, y, Mann J || #4 3)/2 


Under the conditions put on ¢, J*¢ can be shown to be an analytic function 
of a for rea > 1. Further, it can be shown that J*¢ can be continued 
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analytically into the domain rea > 0 and that [%4 = ¢. Also it can be 
shown that L(I*+1¢) = I*6, rea >0. 


For re a > 3, G@-/2 and VG@-)? vanish on the surface of the cone so that 
equation (14) may be written 


n 


. > [ [ (AVGE-V2— GE-D2V4).v;dS8;, rea > 3, 


27T'(a+1) ¢ 
(15) 


i=1“° 
where the S; are the bounding surfaces of the V; which do not lie on the 


I~d — [=+2f 1 








2¢ 
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; \ r 47 
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7 > Pley. 2 
A VS: oe a 
. 
-- Sy a 
Pa 
cil 
ug 
Fic. 1 


surface of the cone. For the problem under consideration, ¢ = 0, z < 0, 
and also ¢ = 0, vu <0. Therefore, it follows from equation (15) that for 
points P(u, y,z) lying forward of the plane w = z,¢ = 0. Hence V may be 
reduced to the volume bounded by the forecone from P(u,y,z) and the 
planes € = 0 and € = ¢ as shown on Fig. 1 

Since ¢ is continuous, ¢ = 0 on Sj. Further, the conormal to Sj is a 
vector lying on Sj. Therefore 


[| (dVGe-D2_ Ge-D2V¢d).v, dS, = 0, 
S? 
and equation (15) may be written 


U—Zz 


ng = PH4p— 5 | dé j Tee V2) -_99(E, ) dy 


(a+ 1) 
0 y 
1 u—z y+r Ge- v2 
shen cromyeernnnecens l [ £,, 0) -| dn, 
' 2a T(a+1) | dé . (7 = al }t=0 ; 
0 vy-?r 


where 
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Since the integrand of the second integral is of class C® for rea > 3, the 
above equation may be written 





“—" y+r 
ae ] y 2 
I%d [*+2f nT (aal) dé [Gene], o9(€. 9) dn— 
_ a yt+r 
C ; ‘ 
Ia T( x 1) oz dé | d(€, 7, 0)[GO-D?),_9 dy. (16) 


0 y-r 

The right-hand side of equation (16) defines an analytic function of « 
for rea > 0, and, under the conditions placed on ¢, the left-hand side is 
also analytic for rea > 0. Since the equality sign holds for rea > 3, it 
must hold for rea > 0 (see Theorem 3, page 89 of reference 6). Therefore, 
equation (16) becomes, for a = 0, 








s=ry—1 f ae [ 2 


ar J * J f(u—€P=(y— 9)? 29} 
u—z y+r d 
_ t P 
l Co | dé | $(E, 2, ) a . 
J{(u—€)?— (y—)?—2?} 





(17) 
0 y= 
Equation (17) is a differential-integral equation since 4(€, 7,0) occurs 
underneath the integral sign. Integration of the last integral by parts 
with respect to » and taking the limit as z > 0 yields 


; Poe r , . , 
lim dé _$(§,9,0) = —}d4(u,y, 0), 
eo 27 OZ . J ay(u—§)*—(y—n)*—27} a 
0 y r 
whence, from equation (17), 
u ytu-—€ ; 
d(u,y,0) = I L2f],.— dé | - HE.0)ey (18) 
e a Vt(u—&)?—(y—n)*} 
0 y—uté 


Substituting the relation (18) into the last integral of equation (17) yields, 
iter some simplification and interchange of the order of integration, 





ie l [ [ | F(E, 0,0) d€dndl 
DU, Y, Zz) Dm | (os — t)2__ (y—n)?—(z—L)?} 


ct 


a7 J J J Vy (U—S) 





€)®—(y—n)?—@+ 0%} 





l [ [ [ f(é. UE C) d&dndZ 
Vl 








5092,26 
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where V is the volume bounded by the forecone from P(u,y,z) and the 
planes € = Cand ¢ = Oand V’ is the volume bounded by the forecone from 
P(u,y, —z) and also the planes € = ( and = 0. For f = 0, equation (19) 
yields the solution to the first-order problem if g = BS;(Bu, y). 

Equation (19) plays an important part in what follows. It is now easy 
to show that the solution to each linear problem is unique. For suppose 
¢, and ¢, are two solutions to one of the linear problems such that 


bb, = $ 40 
for some point u,y,z in R. Then ¢.(u,y,0) = g(u,y) = 0, and 

[d(u, y,2)]. ” a 0. 
Since both functions satisfy the same partial differential equation, ¢ must 
satisfy $uu—Pyy—F22 = 9 for u, y, zin R. 


Therefore, from equation (19), 6 = 0, which is a contradiction, and the 
uniqueness is proved. 


=~ 


In order to illustrate the method and point out certain limitations, a two- 
dimensional problem for which the solution to the exact non-linear equation 
is known will be considered first. 


7. Two-dimensional wedge problem 
Let the upper surface of the wing be given by 


= eH (3), 


NI 


where ¢ is less than 1/8. 

The solution to the first-order problem in terms of the variables u and z 
® = —f(u—z)H(u—z), (20) 

which satisfies the equation ®,,,—@®., = 0 in a region R composed of the 

points of the upper half-plane less the points on the line uw = z. Further, 


the non-homogeneous part of the partial differential equation for the 
second-order problem vanishes in F# so that the solution to the second-order 


is 


problem will not involve the adiabatic index y. 

Although the solution to each linear problem is unique, new solutions can 
be obtained by considering the limit values of the solutions corresponding 
to certain families of wing surfaces. To facilitate the analysis of the situation 
envisaged, the following function is introduced. 

Let D(«,x) denote a function which has the following properties: 

(i) D(a, x) is a continuous function of x and « for all 2, « > 0, 
2 


... OD eD : ; . ; 
(ii) —— and —— are continuous functions of x and a for all x, « > 9, 


Ox Ox* 


except at x = 0, 
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- . eD : 7 D(x, 0- 
(iii) (lim as x > 0 from negative values of x) — = ) = 0, and 
Ox Cx 
. oD : ae CD (x, 0+ 
(lim as x > 0 from positive values of x) = sl = A, both 
Cx Cx 
for « > 0, and where A is a real number, 
., D(a, 0-) ; ie 7 
iv) —* — 0, and there exists a positive number N such that 
a 
2D(a, 07 , , 
~ ) exists for r>WN, 
CL~ 
(v) lim D(a,x) = 2 H(z), 
, . CL(a, x 
vi) lim ) H(x) 
x0 Ca 
and 
= 2D(a, a ‘ - 
(vii) lim - - (x), where 8(2) = 0, x + 0. 
r—>0 Cx 
There exists such a function D, for example, 
D(x, x) = H(x)[(1—A)a*+1+Az]. 
The relations below follow immediately (a is a positive number): 
a a 
*eD 3 
lim dz =a H (x) dx, 
x0 , OX 
| 0 0 
a a 
. fedDeD 1—A? { . 
lim nl - dx — = | H (x) (x) dx, 
x0 , Cl Cx 2 a 
0 0 
l a 
. feD - 
and lim - 1—A = | &(x) dx. 
10 Cox 
0 0 
[f f(x) is a continuous function for 0 < x < a, then 
a a 
| ey» . I 
lim | f(x) —~ dx (1—A)f(0) = | f(a) 52) de. 
1-0 CL” 
0 0 
The expressions on the right in the above relations are similar to the 
notation of Dirac (7). 
, Consider now a wing whose upper surface is given by 


e D(a, 2). 
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For the case of two-dimensional flow, equation (19) may be written 


zZ t+u-z #(u+ 2) ut+e2—€ 
g=a{ far [ fende+ [ at [ fende+ 
0 t z z 
H(u—2) u-2z-¢ uU—z 
+ f a f neoael— foe ae. cy 
0 t 


The solution to the first-order problem is, since f = 0, 


Uu—2Z 


® = —g f D,(fé,0) dé = —B [ HE) aE = —Plu—z), u>0, 
| 


0 
and, since ® = 0 for u < z, simply 
® = —fB(u—z)H(u—z); 
here D, denotes the partial derivative of D with respect to the first variable. 


In the expansion of the independent variable, equation (9), it is sufficient 
to take x, as a function of u alone. Therefore, from equations (11), 





f(€, f) = —2Bary(é — Bay (é) H(é—L)— 
- ; zi | 
— "Be (y+1) H(é—¢) 8(E—Z) +. (22) 
and g(€) = —H(&) H(BE)+ Bx, (€) 8(BE)—8(E) H(BE)BE } 
Substituting the relations (22) into (21) yields 
f Hutz) i(u—2) 2 C4 u—2 
suz)= 4 [FOae+ [ Foa+ f Foa—fac f peyeae— 
0 Z 0 0 } 
Hutz) ute Hu-z) u-z—C 
— | dl | “aatle ) dé — | dt | Bei(€) ag} +-u—2—Br(), 


: t 0 z 


F 4(] = 
where F(Q) = —2Bx( _ MI ae A*) 





Taking into account that 2,(0) = 0 gives 


lim d(u,z) = 4 ( Fig) dé. 


0 


Thus, in order that the potential be continuous on .%, i.e. on the line 


u = z, the following relation must hold: 
4 
v;(u) = cp (tt MEA (23) 


Whence, 
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which, it is to be noted, is the solution that would have been obtained if 
z= u. The solution, to second order, is 





o(u,z) = | —eB(u—z)+e?(u—z) |H (u—z); (25) 
«M4(1+-y)(1+A) am 
x= U — B | : (26) 


From equation (26) it follows that 
Jy tata) 
4 ' 4p 


Substitution of equation (27) into equation (25) and taking into account 


U 





+- oe)]. (27) 


the scale transformation (6) yields 
€?(%— BZ) _ &&MA(1 +y)(1+-A) 
T iw: 4p4 
[,  <MY1+y(1+A] a) 
x Hz [+ 4B —B%|. (28) 


From the above it can be seen that the potential vanishes along lines 
ire «M*(1+-y)(1+A) 
BL 4p3 
Suppose now that « > 0 and that a shock curve passes through the 
leading edge. The Rankine—Hugoniot relations which must hold at a 





+ are)} x 





given by 





A 


(29) 


stationary shock may be expressed as (1) 


; ;' " a 
(i) q, continuous and (ii) san(~—an); 
¢ 
n 


where q, is the tangent velocity vector, qg, is the component of velocity 
normal to the shock, Ag,, is the change in q,, on passing through the shock 
in the direction of the stream, and a is the velocity of sound. Condition 
(ii) holds if g,, and a are both measured directly ahead of the shock and 
holds with the sign changed if both are measured directly behind the shock. 
If vorticity created by the shock is neglected so that a potential exists, then 
i) is equivalent to continuity of the potential at the shock. Thus, for the 
problem under consideration, condition (i) will be satisfied if the shock is 
taken as one of the lines given by equation (29). 

Let @ denote the angle the tangent to the shock curve makes with the 
z-axis. Then, from equation (28) and condition (ii), 


1 «M4(1+y) 
© hana aan cD, 30 
hr” ited, (30) 
From equation (29) it can be seen that the line along which the potential 
vanishes for A 0 is the same as the shock curve determined to first order 
in €. 


tan 6 
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For A = 0, the pressure coefficient is given by 





ge Be, (y- +1) MA 4p? , 
Cp = — a= B+ 2p! o"« 


If the surface of the airfoil is given by 





z= e#+8(#)]H(a), 


oI) 


> 0, 


where s(0) = s’(0) = 0 and s(#) has a Taylor’s series expansion fcr % 
it can be seen that the second term does not affect the determination of z,(u) 
but that it will affect the pressure at the surface of the airfoil. The pressure 
coefficient for this case is 





%e M4— 462 — 
C, = B [1+-8’(x)|+¢? ee #6 [1+<8’(x)}*, (32 


which is the second-order approximation obtained by Busemann by a differ- 
ent method (8). This formula for the pressure coefficient was also obtained 
by Van Dyke (5). 

For the case of « < 0, the solution to the exact equation is the well-known 
Prandtl—Meyer flow of expansion around a corner. The potential in polar 
coordinates for points (p, @) inside the fan is of the form 


$ = OC, psinC,8, 


where C, and C, are constants independent of e. Thus it can be seen that 
the form of the solution assumed, equations (9), is incorrect for the region 
inside the fan which is bounded by the free-stream and downstream 
characteristics which pass through the origin. Using the first-order solution, 
it can be shown that, to first order in e, the downstream characteristic 
passing through the origin is given by 


i= +. (33) 


2B 
The free-stream characteristic passing through the origin is given by 


a (34) 


B 
But equations (33) and (34) are the same as equation (29) for A = 1 and 
A = —1 respectively. That is, the values of A= 1 and A = —1 when 
substituted into equation (29) give the boundaries of the fan to first order ine. 
The values A = 1 and A = —1 are also associated with two families of 
airfoils whose slopes at the leading edge are respectively —e and « and 
where each family approaches the line Z = ef (e < 0, > 0) as a, the 
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parameter, approaches zero. From equation (21) it can be seen that the 
proper value of A for determining the pressure at the surface of the airfoil 
1s zero. 

The above method will now be applied to some three-dimensional wings. 
8, Wing with straight leading edge rolling at constant angular 

velocity 

Consider a flat wing of span 26 and chord c with each tip raked outside the 
Mach cone from the end of the wing. Let the wing be rolling with constant 
angular velocity, w, about an axis parallel to the free-stream velocity vector 
and lying in the plane of symmetry as shown in Fig. 2. 








Fic. 2. 


The appropriate expression for the upper surface of an equivalent wing 


giving the same pressure distribution is 


z = ef9,Z,7 in W (see Fig. 2), 


where € — =, and [el< ‘ 
U ' ~ bB 
The solution to the first-order problem in terms of the w, y, z variables is 
@ -By(u—z)H (u—z), (35) 
which, when substituted into equations (11), yields 
Puu—Pyy — Pez 


— 2 Fh) Bys(u 2)+B|? Sl4y ain) wall H(u—2)— 
ou oy oy" ou* 
M*y*(y- ee =P) + 2(u—2)H(u—2)], (36) 


and =f. (u, y, 0) y?H (u)+ Bx, (u, y)d(u)—B?u? A (u)+ up y8(u). (37) 


— M? 
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Substitution of equations (36) and (37) into equation (19) and taking the 


limit as u + z from values of u > z 2 gives, as a condition for the continuity 
of the potential, that 


att 4. “he 
x,(u,y) = ae ro A) say, (38) 





From equations (9) and (38), 


ss ae  MA(1+y)(1+A) | 


4p3 


For A = 0 the expression for the second-order potential is 








’ [2M2(é—2)+4€K] dédndt 
Hows <5 | tee 
Af E [22 +4€K] dédndt 
JJ. Het —(y—n)?—(2+¢)7} 
Uu—Zz yt+r 





+: f a | a; (n°-+B%€?) dy 


m {(u—€—(y— nF} 
0 y-? 


» (40) 
M(1+y) 
B? 


Integrating and collecting terms gives 


#40, ys2) — [AD ~~ (ee 4 2A) 


8 





where Kk = and r= J{(u—£)?—2?%. 


97 


ae 


- (w—a)(y2— =) H(u—z). (41) 


Using equations (35), (39), and (41) and changing to the Z, 7, 2 variables 
gives for the complete second-order potential 


ME, 9,2) = |- Pe (ee (@ 7 an a" = (e+ | a} 


+222) (pe EPP) Aa (2 aH |. (42) 
Be \7 8 4p T 4p . 


The pressure coefficient at the surface of the airfoil is given by 








C, = — 2,(#, 7, 0)+ O(€?) = + | ap — Pass s+ j(1—2K)]. (43) 
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THREE 


Perhaps a more interesting case is a wing of the same planform whose 
upper and lower surfaces are given by 


i= ea —=),2 ,g in W. (44) 
c 


If this wing is rolling at constant angular velocity w, then the appropriate 
expression for the upper surface of the equivalent wing is 


i c(ta7+ Le — =) (45) 
- 
where ke = —F 


The potential is found to vanish on the surface 


eM4(1+-y)(1-+-A(1+ky) 
7” 4p | 


and the pressure coefficient is found to be 


(46) 











2¢ 2¢\ 22 [x coke 3k*e 2K 
( Fe a oe 
p ( kj : aa ( 5 6p? a i -K) + 
| (4—K)(1+k9) _ 4a 
-- ; +n — ‘ 
Setting Z/c = x, 7/b = y, and —keb = p, the helix angle, gives 


2 1 [a2c2p? 3K 
Oy = 5le—22)—pu — ae | (1 — ge) + 22K 
id 2 apie . (47) 


A plot of C,, against the span is shown on Fig. 3 for M = v2,p = e = 0-1, 
and b/c = 2. 

The rolling moment coefficient due to rolling is given by 

. moment ) ma ep ; - 
C, a F(r+i(e+)—s _ (r—1)(3r2+4r+ 3), 
pune? — 308 
where S is the area of the wing, p is the free-stream density, and r = d/b. 
The ratio of the second- to the first-order term in the above equation is 
R (C,,) end order 0- 4e( (2p°— 3)(r—1)(3r?+-4r+-3)_ 


ace 5B (r+1)(r?-+-1) 
For r 4 and M = v2, 
R= —0-124e, 
which is quite small. The reason for this is that the moments from each 
side of the wing due to the second-order terms tend to cancel out. 
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9. Other planforms 


The method is applicable to other planforms. For example, a triangular 
wing with supersonic leading edges can be treated. The pressure distribu- 
tion to second order can be readily found for the area lying forward of the 
trace of the Mach cone from the vertex, but for the area lying inside the 
Mach cone from the vertex a practical difficulty arises in the evaluation of 
the integrals involved. The type of difficulty encountered is illustrated by 
the following example. Consider a rectangular wing whose upper and lower 
surfaces are given by 


NI 


= +797H(%)H(9). (48) 


Such a wing violates the condition that the surfaces have a Taylor expansion 
for all 7. However, since the expression for each surface is a continuous 
function of 7, the expansion for the potential can be carried out to the 
second-order term, i.e. 


b= «D+ 24+ Ole). 


In terms of the w, y, z coordinates, fory > 0,2 > 0,andu > ,/(y?+2?) =r, 
the first-order potential is 


© = —$Py(u+r—2z)-+ (By/m) | sin{y(g—2")-}} dé +(B/x) [ (2 —r*) de. 
u u (49) 


The above integrals can be readily evaluated in terms of elementary 
functions, but for the purpose intended it is convenient not to do so. From 
equation (49), 


©, = —$By[1+(2/7)sin-{y(w?—z*)-#} ]— (B/z)(u?—1?)!, (50a) 
% = —4p(ut-r—2z) 48 [ sin-rgyie—24)-4 dé, (50b) 
7 


®, = 4By+(By/7)tan—u{|yz—| (u2—r?)-*} — (Bz/7) ( (€2@—r2)-* dé, (50) 


v 





Pum = SCN CE eat (eee 

®,,, = —3B8—(B/7)sin-{y(u?—z*)-4}, (50e) 
and 

o,, = Byte : (50f) 


ue (u®@—2?),(u2—r?) * a/(u2@—1?) 
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For y < 0, u > 1, the first-order potential is 


Whence 

©, = —(By/m)sin“fy(u2@—2*)-}}— (B/m)(u®—r®)!— By, (52a) 

e 

®, 4B(u ry 4° sin-{y(é2—2?)-*} dé, (52 b) 

©. = By+(By/7)tan-lu{|yz-}| (uw? —r?)-*#} — (Bz/7) [ (€2@—r?)-? dé,  (52¢) 
® Bur/(u" 2. (52d) 

me (u?—z?) 

0,,, (B/7)sin-Hy(u* z*)-t_ 48, (52e) 
and 

o,, — Ae. (52 f) 

” ar(u2—z?) 


The second-order potential equation may be written 


M?K, ® CX, (U,Y,z 2M? Ox 
sail ae ®,,,, J - 1 Py 19 CX (UY 2) —®,, meal @,+-2 “1 ams 
? ‘ si B? ou ‘ 
2M*®, xr, x, x 
ae 1 pet 1 = 
—9,, j ee +2 a + | sat St — » (53) 
where K, = y—1+(y+1)B? 

Following the technique of Lighthill (1), the singularities on the Mach 
cone (u2—y?—z* = 0) which are of the form constant x (u?—r?)-? will be 
removed from the non-homogeneous term of equation (53) if the coefficients 
of ®,,, and ®,. are made to vanish on the Mach cone. For y > 0 this can 
be accomplished by setting 


M?K,u M? 
4(U, Y,Z) o( a) o( B ‘), 





where o* = lim 0, = — By, 
u—r 
and o* — lim®. = By. 
i u>r 

f _ Me “Ki, M? 

Whence %, 1 wy — —_ 2y. (54) 
2B B 

For y 1 vy 0, 

since ©, = ©, = ©, = 0 on the Mach cone. 
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Lighthill shows (1) that the determination of the shock strength involve 
solving for x,. This will not be done here. 


For y > 0, u > r, the non-homogeneous term is 


®, M?K, ®.\/(2M*\ 2M? 
fuse) = Onl Gu) “ap + Duly —B) ar) pr Ow Py 
For y < 0, u > r, the non-homogeneous term is 
f(u.y,2) = ~"z (K, ®,,,, P,+2®,,0,+20,,, 0,). 








~ 2 50 75 100 





‘. . r % Semi-span 
~ Sp, 
Fay, - 0-4 Bi Cee, 
~. L Ce, 
M= v2 i S2 
~ 
Helix angle = 0-1 radians >," hag 
by “ ™ Pay. 
4 =2 Ney, 
x Intersection of the chordwise -03 Lm 
station with the end of the wing : ‘ 
ie! 
Fic. 3 


For y > 0, wu < r, the non-homogeneous term and 2, are given in § 8. 
For y < 0,u <7, 
D(u, y,z) = d(u,y,z) = 0 

The expression for the second-order potential is obtained by substituting 
the above values into equations (11). The resulting expression is quite 
complicated and the matter will not be pursued further in this paper. 
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SUMMARY 


It is shown that when an articulated quadrilateral moves in its own plane in such 
a manner that the normals to the paths of two opposite vertices meet on the diagonal 
through the other pair of vertices, then the normals to the paths of the latter vertices 
meet on the diagonal through the first pair of vertices. This is also true for a spherical 
quadrilateral. The properties of the plane and spherical Peaucellier cells are special 
cases of the property. 


1. Introduction 

THIS paper gives a statement and proof of a property of a general plane 
articulated quadrilateral moving in its own plane; the characteristic 
property of the plane Peaucellier cell is a simple special case of this. The 
same property is true for a spherical quadrilateral moving on the spherical 
surface, and the property of the spherical analogue of the Peaucellier cell 


is a particular case of the general property. A generalized interpretation 
of the property is given in section 6. 


2. The basic propositions 

The proof of the kinematic property discussed here depends on the 
geometrical theorem of Pappust and on the following well-known kinematic 
theorem: when three rigid bodies move in the same plane, the three 
instantaneous centres for their relative motions when taken in pairs are 
collinear at any instant. In the sequel we shall identify bodies by numbers 
and shall adopt the convention that rs is the instantaneous centre for the 
bodies r and s. In accordance with the theorem, the three points rs, st, and 
tr are collinear. 


3. Statement and proof of the property 

The quadrilateral, which may be ‘open’ or crossed, consists of four rigid 
and straight links which lie and move in a plane and which are pin-jointed 
to their neighbours; the lengths of the links are arbitrary, subject to the 
quadrilateral closing and its sides possessing freedom of relative movement. 
In accordance with the system explained in section 2 we allot the number 


+ See, for example, H. F. Baker, Principles of Geometry, vol. 1, chap. I. 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 2 (1954)] 
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to the reference plane which is regarded as fixed and in which the quadri- 
lateral moves, while the numbers 1, 2, 3, and 4 are given to the sides (see 
Fig. 1). 

In the figure let the normals to the paths of the vertices 41 and 23 meet 
at P on the diagonal through 12 and 34. Also let the normal to the path 
of the vertex 34 meet the diagonal through 41 and 23 at Q. Then we have 
to show that 12, Q is the normal to the path of 12. 








Fic. 1. Articulated quadrilateral and instantaneous centres. 


Produce 41, 12 and 34, 23 to meet in a point which must be 24 by the 
kinematic theorem of section 2. The instantaneous centre 40 is the meet 
of the normals 41, P and 34, Q, while the instantaneous centre 30 is the 
meet of the normals 23, P and 34,Q. Also the centre 20 is the meet of 
the lines 40, 24 and 30, 23. Now the join of 12 and 20 is the normal to the 
path of 12. Hence our proposition will be proved if we can show that 12, 20 
and Q are collinear. But the points 23, 24, and 34 are collinear and the 
points 40, P, and 41 are collinear. Hence, by the theorem of Pappus, the 
following points are also collinear: 

the meet of the joins of 24 with 41 and of P with 34, i.e. 12; 

the meet of the joins of 34 with 40 and of 41 with 23, ie. Q; 

the meet of the joins of 23 with P and of 40 with 24, i.e. 20. 


4. Extension to the spherical quadrilateral 
Both basic propositions of section 2 are true when we substitute spherical 
surface for plane and great circle of the sphere for straight line. Hence the 
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spherical quadrilateral shares the property of the plane quadrilatera| 
enunciated in section 3, but it must be understood that the word ‘normal’ 
stands for ‘normal great circle’. 


5. The Peaucellier cell and plane and spherical inversion 

Let the quadrilateral become a rhombus ABCD (see Fig. 2) and let the 
opposite vertices B, D be connected by equal links to a point O fixed in the 
reference plane. Then BO, DO are the normals to the paths of B, D respec- 
tively and they are concurrent at O, which lies on the diagonal through the 


as 


IN 


O 


Fic. 2. Peaucellier cell and normals to paths of A and C. 


remaining vertices A, C. Consequently the normals to the paths of A and 
C meet at N on the diagonal BD. By the symmetry of the figure, the 
angles NAC and NCA are equal and it follows that A and C describe 
inverse curves with O as the centre of inversion. This is the well-known 
property of the Peaucellier cell OA BCD. 

When the figure is plane it can easily be proved that 


rr’ = a*—B?, (5.1) 
where Od me. 
OCcC=r', 
OB = OD = a, 
and AB = BC = CD = DA = b. 


For the spherical figure the corresponding relation, which can easily be 
established by spherical trigonometry, is 
, j 
r r a+b a—b 
tan = tan = tan —__ tan—_, 


- ~ - — 


(5.2) 





where the arcs are measured by the angles they subtend at the centre of the 
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sphere. In accordance with the definition of spherical inversiont A and C 
are inverse points with respect to O. When the sphere is projected stereo- 


graphically with O as pole the projections of A and C become inverse points 


ina plane with the projection of O as the centre of inversion. 


6, A generalized interpretation of the property 

The points 12, 23, 34, and 41 may be taken to be the instantaneous centres 
of four bodies instead of the centres of pin joints. Hence the theorem can 
be interpreted as a general property of the relative motions of five coplanar 
rigid planes, namely the reference plane and four others. 


+ See, for example, M’Clelland and Preston, Spherical Trigonometry, part II, chap. xii. 








ELLIPTIC ELASTIC INCLUSION IN AN INFINITE 
ELASTIC PLATE 
By N. JESSIE HARDIMAN (Westfield College, University of London) 
[Received 13 November 1952] 


SUMMARY 

It is shown that constant stresses at infinity in an elastic plate containing an 
elliptic inclusion of different material induce constant stresses in the elliptic inclusion, 
a result similar to the known fact that a uniform two-dimensional electric field will 
induce in an elliptic dielectric cylinder placed in it a uniform field in the dielectric. 
It is interesting that similar results should hold in the electric and elastic problems, 
since the first is essentially a harmonic problem, whilst the second involves biharmonic 
functions. 

The plate is assumed to be in a state of generalized plane stress, and the method 
consists in finding the complex potentials Q(z) and w(z) inside and outside the 
inclusion which give 

(1) appropriate stresses at infinity or in the inclusion, 

(2) continuity of displacement and mean stresses 7m and 7s across the interface 
of are s and normal n. There will, in general, be an interesting discontinuity 
in the peripheral stress ss round the interface. 

Since these solutions were obtainedt my attention has been drawn to some earlier 
results in this field obtained by Donnell (1) by the Airy stress function technique 
in terms of elliptic curvilinear coordinates. A comparison will reveal considerable 
advantage in using the complex potential. Moreover Donnell’s solutions are all for 
the case in which the rigidities of the materials are the same—an unnecessary 
restriction. Further, Donnell’s treatment does not bring out the simple constant 
nature of the stress fields in the elliptic inclusion. This seems to the present writer 


to be one of the most interesting features of the solutions presented here. 
1. The conformal transformation 
THE z-plane is transformed into a ¢-plane by the equation 

af) = ec(€+AC*), O<A<1], (1.1) 
which maps the slit z-plane on the outside of the circle | A! in the 
complex ¢-plane (Fig. 1). Since 


2'(€) = e(1—Al-*) (1.2) 


the critical points S, S’ of the transformation are the ends z + 2cd! of 
the slit in the z-plane corresponding to the points = +A! in the ¢-plane. 
The unit circle |Z 1 corresponds to the elliptic boundary 


x?/a?+-y?/b? = 1, 
if a = c(1+4), b = c(1—A). (1.3) 


+ The results of this paper form part of a thesis for which the Ph.D. degree of the 
University of London was awarded. 


[Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 2 (1954)] 
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Hence the region | of the elliptic boundary in the z-plane maps on the 
annulus A ¢| < 1, and the infinite region 2 outside the ellipse maps on 
the region outside unit circle in the ¢-plane (see Fig. 1). The elastic elliptic 
inclusion has rigidity », and Young’s modulus £,, whilst the surrounding 
infinite plate has elastic constants yp, and #,. The whole plate is assumed 
to be in a state of generalized plane stress, with zero body force. 


= 6 

















z- plane C-plane 
Fic. 1. 

The complex displacement D = u-+iv and the mean stress combinations 
0=a+yy, D = xx—Vy+ 2ixzy are known (2) to be given in terms of 
the complex potentials as 

8uD = cQ(z)—20'(2)—a' (2), (1.4) 
20 = Q'(z)+0'(3), —20 = -0"(2)+8"(2), (1.5) 

where bars denote conjugate complexes, and 
«+1 = 8y/E. (1.6) 

At a point P on the elliptic interface in the z-plane where the normal 
makes an angle ¢ with the a-axis, the mean stress combinations ©’ = nn-+-8s, 
0’ = nn—8s-+ 2ins are known to be related to the cartesian stress combina- 
tions by the equations 

0’ = O, D' = De-*, (1.7) 

{the element ds of the interface at P in the z-plane transforms to the 

element ds’ in the ¢-plane, the normal to the latter making an angle @ with 


the real axis 


dz idse Up. dv i ds’ e9. 2'(Z) ° eilg—8) 
, ‘ ds 
80 that e2te e2id = (8) —- (1.8) 
z'(C)  1—Ag? 


long the interface |Z 1, using (1.2). 
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2. The complex potentials 
For the inner material in the region 1, take 


Q, = Az, w, = Be’, (2.1) 


where A = A,+iA,, B = B,+B, (if the real and imaginary parts of the 
constants are required). These give (using (1.5)) constant stresses in the 


inclusion @, = A, ©, = —B, (2.2) 


and from (1.4) the complex displacement D, is given by 
8u, D, = (x, A—A)z—2B3, (2.3) 
which is clearly a physically satisfactory solution throughout the inclusion. 
For the outer material in the region 2, take 


Q, = Cz+ Fef-, w, = Gz?+He® log (+ Je2l-2, (2.4) 


where C, F, G, and J are complex constants, with C = C,+iC,, ete., and 
H is taken to be real to give zero resultant couple (2) for the stresses round 
the circle at infinity. 

The constants C, G determine the constant state of stress at infinity, 
and the problem is to express the constants A, B, which give the induced 
stresses in the inclusion, in terms of C and G. 

Now 4(mn+ins) = 2(0’+ 0’) = 2(0+ Me-2'4) using (1.7), and so the 


boundary stresses along |¢| = 1 in the region 1 are, using (1.8) and (2.2), 
pre Sy 4(nn,+ins,) = 24,—2B(1—AL2)/(L2—A), (2.5) 
and using (1.5), (1.8), and (2.4), the boundary stresses along |¢| = 1 in the 


region 2 are given by 
4(1im,+ ins.) = 2C,—[F+24(1—AL2)]/(2—A) + 
+ 0[(1+A02)(H—AF)—2(3—AL)(F+-2J)]/(@—A\(1—AL2)®._ (2.6) 


The continuity of these stresses gives a polynomial identity, the vanishing 
of the coefficients of the only powers which occur (2°, £2, £4, 2°) giving rise 
to four consistent equations, only three of which are independent. Simplify- 
ing and separating reals and imaginaries, these equations give 


F, = 2(A,—C,)+-2(B,—-G@,), KR = —2(B,—G,); (2.7) 
J, = —(1—A*)(B,—G,), J, = (1+A?)(B,—G,); (2.8) 
H = 2(1+A*)(A,—C,)+4A(B,—G,). (2.9) 
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The displacements at points on |f| = 1 in the region 1, using (2.3) and 
1.1), are 

8u, D,/e = €-(Ax, A—AA—2B)4-C(n, A—A—2AB) (2.10) 


and in the region 2, using (1.4) and (2.4), are 
8uy D,/¢ = C-"[«,(AC+ F) AC — 2G1/(1—AL2) + 
moke A2)(ky C—C)—An, F+AF —H] (1—Al?)+ 


+ £3[—Aky C+AC + 2422+ F + 25]/(1—AL?). (2.11) 


Again the continuity of the displacements gives a polynomial identity. 
The vanishing of the coefficients of the powers involved ({-, Z, f°) gives 
three consistent equations, only two of which are independent. Simplifying, 
separating reals and imaginaries, introducing the elastic constants 


ee 
and using (2.7), (2.9), and (1.3), gives the equations 
A, A = C,[a(a+6)?—2y(a?+b?)]—2y(a?—b?)G,, (2.13) 
B,A = C\(a?—b?)(1+- 2y—a)+ G,(a+b)?(1+ 2y), (2.14) 
A, A’ = C,aA’+2G, y(a?—6?), (2.15) 
B, A’ = G,(a+b), (2.16) 


A (a+-b)?— 4aby(1+ 2y—2«)/a, A’ = (a+6b)?—4aby/a. (2.17) 


These equations express the constants A,, B,, Az, B, in terms of the 
constants C1, G,, Cy, Ga, which are known if the state of stress at infinity 
is given. 

If the state of stress in the inclusion is given, then the constants A,, B,, 
A,, B, are known and the equations for the constants C,, G,, C,, G, are, 
together with (2.16), 


WC; A, (1+ 2y)+2yB,(a—b)/(a+b), (2.18) 
xG,(a-+-b)? = A,(a?—b*)(a—1—2y)+ B,[a(a+b)?—2y(a?+5?)], (2.19) 
C5 A,—2yB,(a—b)/(a+b). (2.20) 


When the Young’s moduli are the same, « = | and the only elastic 
onstant needed is y. When the rigidities are the same y = 0 and the only 
elastic constant needed is x. This latter is the case dealt with by Donnell (1). 
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3. Discontinuity in the peripheral stress 

Since 7m is continuous across the interface the discontinuity 88, —&8, 
along || = 1, where £ = e*?, so that @ is the eccentric angle of the point 
of the elliptic interface, is given by 0,—0,. Hence 

33, —38, = [(A,—C,)(1—2?) + 2(B, — G,)(cos 20—A) — 
— 2(B,— G,)sin 20]/(1+-A?— 2A cos 26), 

which can be expressed in terms of the perpendicular p from the centre 
of the inclusion upon the tangent, as 


~ a+b /(F, p? 
88,— 88, = d pits (B,—G,)— 


(B,—G,) 
a—b\2ab 7 


ab 


9 9 


[(a?— p?)( p? b»)]), 





where p® = c?(1—A?)?/(1+A?— 2A cos 26). (3,2) 


4. Plate under given stresses at infinity 
(a) A simple tension 7’ at an angle f to the x-axis at infinity gives, for 
the constants in (2.4), 


C, = f, C, = 0, G, = —T' cos 28, G, = Tsin2f8. (4.1) 
(6) An all-round tension 7 at infinity gives 
C, = Sf, C, = G, = G, = 0. (4.2) 


(c) A simple shear S at infinity gives 


C= ¢ = @, G, = 2Ssin 24, G, = 28 cos 2¢. (4.3) 


5. Given state of stress in the inclusion 
(d) A simple tension 7 at an angle f to the x-axis in the inclusion gives, 
for the constants in (2.1), 
a, oF, A, = 0, B, = —T cos 28, B, = T sin 2p. (5.1) 
(e) An all-round tension 7 in the inclusion gives 
A, = 27, A, = B, = B, = 0. (5.2) 
The special cases in sections 4 and 5 are discussed in greater detail in the 
writer’s Ph.D. thesis (3). 
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MOTION OF AN ELLIPSOID IN A 
ROTATING FLUID 


ON THE FREE 





By K. STEWARTSON 
(Department of Mathematics, The University, Bristol) 


[Received 27 November 1952] 


SUMMARY 
It is shown theoretically that, if a spheroid of density po is placed symmetrically 
n the axis of a rotating fluid of density p, it will be in stable equilibrium for 


translational disturbances if o 1, and for rotational disturbances if o > 1 in the 
ease when it is oblate, or if a 1 in the case when it is prolate. Further, if 1 > o > 3, 
t will have a single period of translational oscillation, while if o < } it will have 
two such free rotational periods ; and if stable, it may have also a free period of rota- 


tional oscillation. Some experiments on the free motion of bodies in rotating fluids 
wre also described. These are in substantial agreement with the theoretical predic- 
tions. It was also observed that there are asymmetric positions of equilibrium, stable 
with respect to angular displacements, having the axes of the bodies horizontal. 


1. Introduction 

[x previous papers (1, 2) an account was given of the uniform slow motion 
of an ellipsoid in a rotating fluid, when one of its principal axes was paralle] 
to the axis of rotation. The ultimate motion of the fluid relative to axes 
rotating about the axis of rotation with the same angular velocity as the 
undisturbed fluid is steady but nevertheless has some remarkable properties. 
In a description of some of them which is given below it will always be 
understood that we are discussing the velocity of the fluid relative to these 
rotating axes. 

Firstly, the motion is entirely two-dimensional, the components of the 
fluid velocity at any point depending only on distance from the axis of 
rotation and on the azimuth angle with respect to the same axis. Secondly, 
the cylinder C circumscribing the ellipsoid and having its generators parallel 
to the axis of rotation separates out two regions with markedly different 
characteristics. These characteristics may best be described by considering 
the motion of the ellipsoid along and at right angles to the axis of rotation 
separately. If the ellipsoid is moving along the axis of rotation, then the 
fluid inside C'is pushed along in front with the same velocity as the ellipsoid 
and, in addition, the streamlines relative to the ellipsoid are ellipses: outside 
( the fluid velocity is parallel to the axis of rotation, is singular on C, and 
tends to zero as the distance from C tends to infinity. If the ellipsoid is 
moving at right angles to the axis of rotation the streamlines inside C 


(Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 2 (1954)] 
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relative to the ellipsoid are plane sections of congruent ellipsoids, but 
outside C the relative streamlines lie in planes at right angles to the axis 
of rotation: further, in those planes the streamlines have to execute sharp 
turns either as they near C or just after leaving it. 

Xemarkable as these results are, it is even more remarkable that the 
agreement between theory and experiment should be so good. Both for 
motion along and perpendicular to the axis of rotation the cylinder C was 
observed (3) to separate out two regions with markedly differing charac- 
teristics and in the first case a cylinder of fluid was seen being pushed in 
front of the sphere. In the second case it was observed that outside C the 
streamlines on one side broke away from C just before reaching it, and 
that the fluid inside C had no motion relative to the sphere, but both of 
these discrepancies can be explained by the known reluctance of real fluids 
to execute sharp turns. 

Nevertheless the ultimate behaviour of the fluid is most unusual, and it 
is of interest to inquire both theoretically and experimentally whether the 
fluid would behave in this peculiar way if the body were free to move, under 
gravity, for instance, and not, as in the previous work, constrained to move 
with uniform relative velocity. We ask in fact whether such a motion could 
be stable. The motion of translation of an ellipsoid in a rotating fluid was 
discussed in (2). It may be deduced from the results obtained there that, 
if U; and V; are respectively the components of the velocities of the body 
and of any point of the fluid relative to a set of axes rotating with the fluid, 
then 


ytot 
: l i 
ee | 
2m | 
y—ta 


U,T,;e"ds, y > 0, 


where U, = | e-*U; dt, 


and 7;; is independent of U; and regular, apart from branch points with 
exponents —3 on a strip of the imaginary axis of s. Thus, any instability 
of the fluid in translation can only arise from poles or branch points of U; 
with positive real parts, and hence is due to an instability in the motion of 
the body itself. A similar result may be shown to hold when the body is 
rotating relative to the fluid. It is therefore sufficient to consider the 
stability of the motion of the spheroid itself, and further, since all relative 
motions are assumed to be small, we may, in such an investigation, neglect 
all the external forces. The effect of an external force such as gravity is 
merely to superpose a steady motion on to the free motion. 
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We consider in particular the stability of a spheroid with principal axes 
of lengths 2a, 2a, and 2c, and whose axis of symmetry is parallel to the axis 
of rotation. The spheroid is supposed to be in relative equilibrium in this 
position and, owing to the action of viscosity, rotating about the axis of 
symmetry with the same angular velocity as the fluid. It is then given a 
small displacement and allowed to move under the action of the pressures 
on its surface only. The pressure distribution is determined as in (2), and 
on writing down the force and couple on the spheroid and equating them 
to the rates of change of linear and angular momentum respectively, two 
transcendental equations are obtained. The motion of the spheroid is stable 
only if the roots of these two equations are all real. 


2. Results 

The force exerted by the fluid on the spheroid may be divided into two 
parts. First of all the fact that the spheroid is not symmetrically placed 
with respect to the axis of rotation means that there will be a force exerted 
on it, tending to move it towards this axis, and which we may call the 
centrifugal force. Secondly the motion of the body relative to the fluid will 
induce a force on it, and this we may call the hydrodynamical force. 

We find that if the spheroid is displaced bodily along the axis of rotation, 
the only force acting on it is hydrodynamical, tending to bring it to rest. 
If the spheroid is displaced bodily in a direction at right angles to the axis 
of rotation, the stability is dependent solely on o, the ratio of the densities 
of the spheroid and the fluid. If o > 1 the spheroid is unstable, if o < 1 
the spheroid is stable and oscillates with one or more free periods of oscilla- 
tion, while if o = 1 the spheroid ultimately comes to rest relative to the 
fluid. The instability when o > 1 is due entirely to the hydrodynamical 
forces: if the centrifugal forces alone were of significance the spheroid would 
be stable for all values of c. The special case when o = 1 and a = c has 
also been considered by Grace (4, 5) and we are able to confirm his results 
as special cases of ours. 

When the spheroid is given a small angular displacement about its centre 
we find that the resulting motion is unstable if (s—1)(c2—a?) > 0. Other- 
wise it is stable and either oscillates with a single free period or comes to 
rest relative to the fluid. Ifo = 1 the spheroid also comes to rest ultimately, 
while if c = a a rotation of the sphere does not affect the fluid. Again the 
hydrodynamical forces play a large part in determining the criterion, for 
ifthe centrifugal forces only were of importance the condition for instability 
would be o c4/a4—1. 

We conclude, therefore, from theory, that the only stable spheroids with 
their axes parallel to the axis of rotation are prolate ones lighter than the 
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fluid they displace. Prolate spheroids heavier than the fluid they displace 
are unstable with respect to translational and rotational displacements, 
Oblate spheroids heavier than the fluid they displace are stable with respect 
to rotational displacements but unstable with respect to translational 
displacements. Finally, oblate spheroids lighter than the fluid they displace 
are stable with respect to translational displacements but unstable with 
respect to rotational displacements. 

These conclusions are substantially in agreement with the experimental 
observations described in the last section of the paper. It was also 
observed that if c—1 were just negative (say about —0-01) there was a 
slow steady motion parallel to the axis of rotation even if the symmetrical 
position were unstable, for there was in addition an asymmetric position of 
stability with the axis of the body horizontal, and therefore perpendicular 
to the axis of rotation. This motion may be explained by means of the 
results of (1), for gravity will exert a resultant force upwards on the body 
so that it would move upwards with constant acceleration were it not for 
the hydrodynamical resistance which in a slow motion is ultimately propor- 
tional to the velocity. A state of equilibrium will clearly be reached in which 
the body moves parallel to the axis of rotation with uniform velocity. When 
o—1 was just positive a similar motion was observed except that this time 
the bodies moved downwards and the translational instability gradually 
forced them near the side of the beaker containing the fluid. Further, when 
o > 1 there was always an asymmetric position of stability with the axis 
of the body horizontal. When there were two positions of stability with 
respect to angular displacement, as in the case of oblate bodies with o > 1, 
they preferred the asymmetric position, in which their angular momenta 
were least. 

Let it be said, however, that the experiments were not completely satis- 
factory because the angular velocity of the fluid were so small that o had 
to be kept very near to unity, which meant that instability was not marked, 
and also that any non-uniformity in the bodies assumed an importance out 
of all proportion to its size. It is hoped therefore to repeat the experiments 
with a more rapidly rotating fluid. 


3. The hydrodynamical forces 

We suppose that the fluid is rotating about the axis Oz with angular 
velocity 4, that O is a fixed point on that axis and that O(x,y,z) form a 
right-handed orthogonal system of axes in which Ox, Oy are rotating about 
the z-axis with the same angular velocity as the fluid. It is supposed that 
the spheroid is initially at rest relative to these axes, that its centre is 0, 


and that c2(a2+ y?)+ a2? = are? (3.1) 
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is the equation to its surface. We restrict our attention to a symmetrically 
disposed spheroid because the stability of such a body can be fully investi- 
gated; a similar formal solution may be found for a general ellipsoid, but 
itismuch more difficult to obtain useful information from the transcendental 
equations. 

In the motion, let the spheroid undergo a small displacement in which 
the centre moves to €, », ¢ and the axis of symmetry of the spheroid has 
direction cosines (1, m, 1) relative to Oxyz, where €, y, ¢, 1, and m are all 
small, and let it move freely after a given initial motion. The equation of 
the spheroid is then in general 

F = c7{(a—£)*+(y—n)?+(z—C)?]+ (a? —c?)[z—C+lx+ my}?—a*c? = 0. 

(3.2) 
Ifu, v, w are the components of the fluid velocity relative to O(x, y,z) and 


are all small, then the boundary condition D¥/Dt = 0 on the spheroid 


reduces to 


9 9 
Ux vy. W2 Ze y.,2, @&—ce . : a 
L | €1+—9+—C ——.— (ral+- yzm) (3.3) 
a a“ c* a* a* c* a“c~ 
> 9 9 9.9 9.9 
when c*(2*-+-y*)+a*2* = a*c’. 


The boundary condition at infinity is that uw, v, and w should tend to zero 
there. Since wu, v, and w are all small, at least to begin with, we may 
linearize the equations of motion by following previous writers on rotating 
fluids, notably Morgan (6). It is important to bear in mind, however, that 
such a linearization, while justified in practice, cannot be defended on 
mathematical grounds when ¢ is large, because the vorticity, and hence the 
terms neglected in the linearized equations, oscillate infinitely. The agree- 
ment with experiment is sufficiently good, however, for us to assume that 
no serious errors arise from this cause. 

The boundary condition when t = 0 is complicated by the impulsive 
motion of the spheroid at that time. One plan is to follow Morgan (6, 7) 
and assume that at ¢ = 0+ the relative motion of the fluid is the same as 
ifthe fluid were at rest at infinity. Morgan showed that the only difference 
from the solution when the fluid is assumed to be at rest at infinity relative 
to the rotating axis is that the pressure term no longer contains the delta 
function describing the impulsive pressure in the fluid. The disadvantage 
of Morgan’s method is that the boundary condition on the spheroid is 
complicated. Instead we shall suppose that at ¢ = 0— the fluid is at 
relative rest and at that time an external impulse is applied to set the 
body in relative motion. Hence we include in the pressure terms a delta 
function corresponding to the impulsive pressure in the fluid at t = 0: no 
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complication arises from this, and we obtain the same solution for t > 0 as 
we would with Morgan’s initial condition. 
If p is the pressure and p the density, we may write 


P = i“ L(a2+y?), (3.4) 


Further denote by a bar the Laplace transform of any function with respect 
to t, so that, for example, 


0 ytio 
- , a 
P= |e*“Pdt and P=— e"Pds, y>O0. (3.5) 
21 
0 y — too 
The equations of small relative motion are (6) 
Cu oP Ov. oP 
—v=- — ; . — + = — — 
ct Ox ot oy 
Cw oP ou ov, ew 7 
=——, and —+—4—= 90. (3.6) 


ct oz Cx Oy Oz 


Hence we have, applying the Laplace transformation with respect to f, 





~~ eP aP ee eP aP - eP 
(1+s*}8 = —s— — —_, (1+s*?)6 = —s—+—, and s#=-——, 
ex = ey ey ex ez 
(3.7) 
2D 2D 2o2)D 
O° oF 1-+-s? 67} 
so that —+—++—__ —_ = 0. (3.8) 
ox* Oy?” a @f 
In terms of P the boundary condition on the spheroid is 
2é@P yeéP 1+8%2zéP xeP yéP 
Be a — ye See emma. fl aafien demi Sei. smal 
a* cx a* Cy 8? c? dz] at dy a? Ox 
L,s I oa , 278 
a (8s§—§o) + a (8}—) + a (s¢—Lo)— 
(1-++s?) » (3.9) 


2__ 2 
— — [az(sl—1,)+- yz(sm—my)] 
where the suffix 0 refers to initial values; further, since uw, v, w—> 0 as 
x, Y, 2 > ©, it follows that grad P > 0 as well. 

If (X,Y, Z), (A, B,C) are the components of the resultant force and 
couple at the centre of the spheroid exerted on it by the fluid, then 


X+i¥ = —1M(E+i9)—p [| P(f+ig) as, 


es) M (c?—a?)(l+im)+ip [ i Plh(xa+iy)—2(f+ig)] dS, 


20 


Z=p|{PhdS, C=0 (3.10) 


; 
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where f, g, are the direction cosines of the outward drawn normal to 
and the integrals are taken over the surface of the undisturbed 
spheroid, and M is the mass of fluid displaced by it. Let the external 
impulse required to start the motion have components (Xo, Yo, Z)) and a 
moment about the origin with components (A,, By, 0). Then the Laplace 
transforms of the applied externa! force and couple have components 
(Xo, Yo, Z») and (Ay, By, 0) independent of s. The simplest way to calculate 
these is to note that they are just sufficient to counteract the impulsive 
pressure of the fluid at ¢ = 0, and therefore 

lim (X+ X,) = 0, ete. (3.11) 


8% 


The general nature of the solution for P has been discussed in (2) and 
here we need only write down the appropriate form. It may be noted in 
passing that the branch points of the solution all occur on the imaginary 
axis of s in the strip —i < s < iso that the stability of the fluid is linked 
with the stability of the body unless the spheroid oscillates with a period 
of greater than half the period of rotation of the fluid, when resonance 
would oceur. It turns out, however, that the periods of free oscillation of 
the spheroid are all less than this. 

It is convenient to consider the translational and rotational disturbances 
separately. It will become clear in the course of the argument that they 
are independent. 


4. Translation (/ = m = 0) 

The differential equation satisfied by P may be reduced to Laplace’s 
equation with the same spheroidal boundary if we replace z by zs(1+-s*)-* 
and ¢ by cs(1-+-s?)-? so that the problem is similar to that of the translation 
of a spheroid in a fluid at rest at infinity. The only difference, in fact, is 
that the boundary conditions are now slightly more complicated, but 
the same general approach may be used. The details of the method have 
been given in (2), and here we need only sketch the steps taken and 
write down the results. We introduce spheroidal coordinates (A, », @) 
orthogonal with respect to (3.8), of which A is zero on the spheroid. Then 
we find that A satisfies 
— 2 
+ | @-2-wy-Fea|=% 

| g2 l+s? 1-+8? 
tends to infinity with s, and qua function of s has branch points on the 
imaginary axis of s between +7. 

From the boundary conditions on P it is clear that we must look for 
solutions of (3.8) which are of the form 2, y, or z multiplied by a function 
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of A and which tend to zero as 2?+-y?+z? > o. The appropriate solution 
is then 
7 dr , dr 
x + By) ToT ET Re de | ree, §€=6 
Jr +IPRe te” J rt IG+) 


Naty Nap 
in which «, 8, and y are constants and % = c?s?/a?(1+-s?). Substituting into 
the boundary condition on the spheroid, on which A = 0, and writing 








H, = | ee o_o eT 

J (br +1)(7+1)! © Jd (or +1)(7+1)! 

0 0 
we get, on substituting in the boundary condition on the spheroid, 

x(sH, —2s)+ BH, (s?+-] )(s& —E,) 

and — v1, - B (sH, — 2s) — (gs? + 1)(s7— 7), (4.4) 
whence r+i8 — (ott = — a fo— to] (4.5) 
Further y = ee 5 (4.6) 


Thus we have expressed the pressure in terms of the displacement of the 
spheroid. We can now determine the motion of the spheroid by finding the 
force which the pressure exerts on the spheroid and equating it to the 
acceleration of the body, thus obtaining three equations for €, 7, and £. 
The components of the force exerted on the centre of the spheroid are 
(X,Y,Z), where 
22 
X+s¥ = —pmtiqn STV foes iq)—e,—in, (4.2) 
2s—(s— 7 


= Ss MH,(s¢ -Ca) 
and Z = — 20) 4.8 
An 5 Tl, (4.8) 





while the couple vanishes. 

The translation of the spheroid along and at right angles to the axis of 
rotation has therefore no effect on the relative rotation of the spheroid. 
We shall see below that the converse is also true (5.4). 

We consider the motion along the axis first. Let the density of the 
spheroid be po so that its mass isoM. The equation of motion is independent 
of € and 7» and is Mot = Z+%, (4.9) 


where Z is defined above and Z, is the impulse required to set the body in 
motion. Since f = f), = £, when t = 0, 


_ MsH,(sf— lo) 


Mo(s?{—sl,—£,) - 4-Ey (4.10) 
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Now Z, is an impulse function at ¢ = 0 and therefore Z, is independent of s. 
Hence letting s + 00 we find that 
Zo MC, 


Hs | (4.11) 


H,—2 


and is the same as that required to set the spheroid in motion if the fluid 
were at rest at infinity, in agreement with Morgan’s result (7). The displace- 


ment of the spheroid is 


1" Zy+Mot,)(H,—2) e 5. 
271M (o—1)H,—2c_ 3? 


(4.12) 


on inverting the transformation. 

It may be shown, using methods similar to those in the Appendix, that 
;>+H,(2—H,)-* has no zeros in the s-plane cut along the imaginary axis 
from —i to 7. Hence the ultimate value of € may be found by examining 
the pole of € at s = 0 and the nature of the branch points at s = +i. We 


: , , Ws Hi, cost ; 
¢ = b+: blots ay +9 >) (4.13) 


find that 


when ¢ is large, and so when ¢t = 4z, i.e. after one revolution of the fluid, 
the spheroid has very nearly come to rest. In particular, when o = 1 and 


c the sphere comes to rest at 
C=%+—b, (4.14) 


in agreement with Grace’s result (4). 

Hence, we concludethat not only isthe spheroid stable under disturbances 
long the axis but that the restoring forces are dissipative and ultimately 
bring it to rest. The motion described in (5) can therefore be attained by 
allowing the spheroid to move along the (vertical) axis of rotation under 
gravity provided only that the fluid is rotating fast enough. This has been 
confirmed experimentally and may indeed be observed when stirring tea. 
So long as the tea is rotating rapidly the tea leaves appear to remain almost 
stationary relative to the fluid and it is only when the rotation has nearly 
ceased that the tea leaves fall to the bottom of the cup. 

Let us now consider the motion at right angles to the axis of rotation. 


Combining the equations into one we obtain 


Mo(€+ij—H+ié—}é—}y) = X+iY4+X,4+tY). (4.15) 
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Hence since € = &,, etc., when ¢ = 0 


Fe acca alee IEh | 
i to t+18)- 4 (s—i)H,— 2s] 





(€+47) 
s(s?+-1)H, 
ie ti aw 
(o—i)H, — | —— 
As in the case of Z, above, X, and Y, are independent of s and hence, 
letting s > 00, we obtain 
v ‘o> : . H 
X,+tY, = M(é,+%,)| —- : 4.17) 
oT *4o (fo is = (4.17) 
The displacement of the centre of the sphere may now be obtained as a 
Bromwich integral, as in (4.12), and it then follows that the motion is stable 
if all the poles of €+-i7 which are the zeros of the coefficient of €+i7 in 
(4.16) lie on the imaginary axis of s. If any occur elsewhere, then by 
symmetry half of them will have positive real parts contributing exponen- 
tially increasing terms to +i and so making the motion unstable. 
We are led therefore to consider the zeros of 
s(s?+-1)H, 
(s—t)H,— 2s 





= Xot+iYy+ Mo(é,+ io) + M (Eo+ ing) (e+ t)— 


A(s) = o(s?+is)—}(o—1)— (4.18) 


First of all we prove in the appendix that 
(s?+-1)sH, 
(s—t)H,—2s 
has no zeros or poles in D, a contour consisting of the infinite circle, both 
sides of a cut from s = —i to s = +7 along the imaginary axis of s, and 


small circles round s = 0ands = +7. A similar argument then shows that 
since 








A(t) = —}(9o—1), A(—t) = —}(o—1), and A(i«o) <0, (4.19) 
A(s) has one zero in D if § < o < 1 and for other positive values of o has 
two zeros in D. Nowif s = ié*, where &* is real and |€*| > 1, —o(s*+is) 


increases with €* and sH,(s?+-1)[(s—i)H,—2s]-! is always positive. Hence, 
if c > 1, A has no zeros, if } < o < 1, A has one zero, and if o < }, A has 
two zeros on the imaginary axis of s. Since the total number of zeros of A 
is known and given above, we may conclude that if ¢ > 1 the spheroid is 
unstable because one of the zeros of A must have a positive real part, while 
if ¢ <1 the motion is stable, with a single free period of oscillation if 


j < o < 1 and two free periods if « < }. If we neglect the forces due to 


9 
the relative motion of the fluid, the criterion for stability may be found 
from (4.18) by neglecting the term involving H,, and is accordingly o > 0. 
Hence the hydrodynamical forces, as distinct from the centrifugal forces, 


are of fundamental importance in deciding the stability of the spheroid. 
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When o = 1, A degenerates somewhat, but we can show by a similar 
argument to that when o 1, that it has no zeros in D. Hence the spheroid 


will ultimately come to rest, and, in fact, when ¢ is large 


, F ‘ ee 2 -  7C) cost 
on v7) Eo r 19 (So io s—77| (—1+Z)+0/ +} (4.20) 
ai l/s=< 


so that the relative motion has almost ceased after two or three revolutions 
of the fluid. In particular when £) = ny = 1%) = 0 and a = c the sphere 


comes to rest at 9 
: OTT ; at 
$ Ss fo; 7 3£o, (4.21) 


which is the result obtained by Grace (5). 


5. Rotation (€ = 7 = ¢ = 0) 

The method for determining P in this case is the same as that sketched 
out in the previous section, except that in virtue of the boundary condition 
in (3.9) we need to search for solutions of (3.8) of the form xy or xz multiplied 
by a function of A. We find that 


r dr 


r Faz+ Gyz ES 
: y) (Jr+1)?(7+1)! 


Alar 


(5.1) 


where F', G are constants and A, s are defined near the beginning of section 4, 
is of the required form. The condition that grad P + 0 as a?+-y?+2? > a 
and the condition on the spheroid are both satisfied, the first automatically, 
and the second if 


(F+iG)|{s+(s—7)}H,— 2s] a b(s?4 1)[s(/+-im)—1,—img], 
c2 
(5.2) 
j dr 
here H ; - 5.3 
where ; | G1) (5.3) 


Further, the moment of the pressures on the surface of the spheroid about 

its centre has components (A, B, 0), where 
= ” m evst — iM e ° 7 Y » 

A+iB 5g he a*)(l+-im) — (c?—a*)(F+iG)H,, (5.4) 

while the resultant force at the centre is zero. This completes the proof of 

the earlier statement that the rotational and translational motions are 

independent. 

The equations governing the rotation are: 


— R 














K. STEWARTSON 


Me (a?+-0*)(l—a— Hl) 4 ato +H) = B+ B, 
v0 vo 
» (5.5) 
and me (a?+-c?)(m+1—1m) — Me a(l—1m) = —A—A, 
5 vo 


where (Aj, By,0) are the components of the impulsive force required to 
start the motion and so are delta functions with Laplace transforms 
independent of s. 

The equations of motion reduce, on taking the Laplace transform with 
respect to ¢, to 
(i+im)T(s) = —i[A,+iB,]+ 


+o(a?+c?)(i,+ ig) + (y+ im)| _ << Ne, (5.6) 





where 


CHa? (c2—a?)? ws(s?+-1)H, 








4 @ [8s +-:h(s—i)]Hg—2s° — 

We may now obtain 4,+7B, by letting s > 00 in (5.6), whence 
-— oe oo i(c?—a?)H, ze 
A,+iB, = — (loin) | ag a 3 (5.8) 


As was the case with +n, the condition that /]+-im should remain small 
is that the poles of /+-im should all be on the imaginary axis of s. These 
are clearly given by the zeros of I'(s) which we now consider, using the 
same method as for A(s). Using a similar method to that described in tl. 
appendix we can show that 











ws(s (s?+-1)H5 _ 
[s+o(s—i) )4,— 
has no zeros or poles inside D. We then deduce that since 
P(0) = —(o—1)(e@—a?) 
if (i) c >a, o >1 ore <a, o <1, T has two zeros in D; if (ii) ¢ >4, 
c?—a? . ; ~—— 
> > => OOree < a,o > 1,T has no zeros in D; and if (iii) ¢ > 4, 
3(3c?+-a?) 
c*—a? ? * * 
I has one zero in D. Further, if s = igé* where €* is real 
~ 3 3c? a?) 2)’ 
and |€*| > 1, —(s+47){(a*+-c?)s+ 4i(a?—c?)} is positive and increases with 
* ¢ 
&*, and ps(s?-+1)H, 
[s+y(s—i)]H,— 


is also positive. Hence, if I'(i) < 0, then I( nee é* > 1, is negative and, 
if (2) > 1, [ has an odd number of roots in i < s < io. Now 


P(t) = —}{30(3c?+-a?)—c?+a?} (5.9) 
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and so in (i) and (ii) above [ has no roots ini < s < io0, while in (iii) it has 


at least one and therefore, since it cannot have more than one, it has one 


7 so, since 
only. Also, si (a2—c2) ; 
C(s) = ———_—. (5.10) 
3a” log(s+-t) 
near 8 i, T has no zeros in —i > s > —ioo. Thus, in (i) the two zeros 


of I are not on the imaginary axis of s and one of them must therefore have 
a positive real part so that the spheroid is unstable. In (ii) and (iii) the 
spheroid is stable, and further in (ii) it will ultimately come to rest with its 
axis along the axis of rotation, while in (iii) it will ultimately oscillate with 
one free period. When a/c = 0, the zero of T may easily be found to be at 


i(1—a) 2 \} - 
‘= satel + (=a) | aes 


If c > } this zero is on the cut and is therefore isolated in that it exists 
when a/c = 0, but when a/c is small T never vanishes. The maximum value 
of the zero of T occurs when o = 0, a = 0 and is 1-2077. 

If c = 1 we need a special argument because ['(0) = 0, but it is on the 
same lines as previously and we find that I never vanishes inside D; when 
o > 1 the zeros of TI tend to and coalesce at s = 0. The spheroid is there- 
fore stable as we should expect and if displaced comes to rest at 


; (37a-+ic)2ac 2c?H,—a?—c? 
l ee LL l 
8m 





+-4m,+ (ly, +7m,)——— 
) ) 0 , er Ere ee r os 
ai 4a°— 8ac? + 37ic*| (a? +-c?)H,— 2a? 


0 
6. Experiments 

A circular glass bowl, 8 inches in diameter and 8 inches high, was filled 
to within 2 inches of the top with water and waxed to a horizontal gramo- 
phone turntable. The angular velocity of the turntable was adjustable up 
to a maximum of 10 radians per sec., but since at higher speeds the free 
surface of the water took a long time to settle down, the angular velocity 
was usually maintained at about 6 radians per sec. 

Taylor’s experiments (3) suggested that a slow motion theory is valid 
only if V <aQ, where V is the velocity of the body relative to the fluid, 
ais a representative length on the body, and Q is the angular velocity of 
the fluid. Any body whose specific gravity is not unity will move vertically 
relative to the water, due to the action of gravity, and using the resistance 
formula obtained earlier (1) we may deduce that the terminal vertical 
velocity will be O(g|1—o\Q-1). Hence, if our theory is to be comparable 
with our experiments, g\1—o| <a?. 
Taking a = 0-05 ft., g = 32 ft./sec.2, and Q = 6 radians per sec., this 
condition becomes |1—o 0-06. The most satisfactory substance with 
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a specific gravity in this range is pure camphor, which has a specific gravity 
of 0-99. Bodies of specific gravity slightly in excess of unity could then be 
obtained by smearing the surface of the camphor with plasticine. Although 
the bodies used were not spheroids, being in fact cuboids, it was felt that 
they were sufficiently close approximations to provide a satisfactory test 
of the theoretical predictions of this paper. 

The vertical motion of the bodies was first examined. It was found that 
if o was near unity the bodies were always stable for vertical displacements 
and when left alone ascended or descended with a strikingly slow velocity, 
When the fluid was at rest the whole vertical motion of the camphor would 
take a few seconds only but when the fluid was rotating at 6 radians per sec. 
the corresponding time was of the order of 5 minutes. An attempt was made 
to compare the calculated and experimental values of the terminal velocities, 
and although the agreement was good, the possible sources of error were 
sufficient to compel us to regard it as being qualitative rather than quantita- 
tive. For example, it was impossible to eliminate entirely the secondary 
motions of the fluid to which the bodies were very susceptible. 

Secondly, the motion of the bodies in a plane at right angles to the axis 
of rotation was examined. Those bodies which were lighter than the water 
tended to be at least an inch from the side of the bow] and most often were 
quite close to the axis of rotation. On reaching the free surface they main- 
tained a position of relative equilibrium some distance from the side. Those 
bodies which were heavier than the water ultimately moved out towards 
the side although occasionally this motion took some time to start; further, 
they were often a long time reaching the side but when in relative equili- 
brium at the bottom of the bowl they were always in contact with the side. 
We may explain this reluctance of the heavier bodies to touch the side as 
follows. If the fluid is at rest then a wall tends to repel a body moving per- 
pendicularly to it with a force that is large when the body is near the wall 
but which rapidly dies away as the body moves away. There seems to be no 
reason why a similar state of affairs should not occur when the fluid is 
rotating, and so in our experiments the side is likely to reduce the instability 
of the heavier bodies. 

Thirdly, the angular stability of the bodies was examined. This was 
more difficult because the bodies had to be perfectly symmetrical, and also 
because there were asymmetrical positions of equilibrium as well, so that 
care was needed to set up the symmetrical position. For example, oblate 
bodies heavier than water were stable when their axes were horizontal and 
it was difficult to get the axis vertical from that position. Presumably the 
reason was that a greater angular momentum was required for a given 
angular velocity when the axis was vertical, so that to compensate for this 
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vity the angular velocity of the body decreased when the axis was made vertical, 
n be and that reduced the stability. If when the axis was made vertical the 
ugh body was given a compensating angular velocity, then the oblate body was 
that quite stable in the symmetrical position. Prolate bodies with o > 1 were 


test } not stable when their axes were vertical and they had an asymmetric 
position of stability with their axes horizontal. Those bodies lighter than 





that the water were found to have a symmetrical position of equilibrium only 
eat if they were prolate. Prolate bodies with o < 1 may have had a position 
ity. of stability with horizontal axes, but it was not very pronounced, while 
vuld the oblate bodies were quite stable in that position. 

ae The stability of the symmetrical positions of equilibrium of the bodies 
ade | — jighter and heavier than water seemed to be lost near the free surface and 


mm, the side respectively, presumably because of the repulsive force from the 


bain bounding surface. 
ita 
ary APPENDIX 
We prove here that 
s?+ 1)sH. 
LXis I (s—i)H,—2s and J A Lada 


(s—i )H, 2s 


te ; i4 re , , 
Wer have no zeros or poles in the contour D, comprising the infinite circle, both sides of 
: I 
ere § the cut from —? to ¢ along the imaginary axis of s, and small circles round s = 0 
in and s 
nia From equation (4.3) we obtain 
x 
rds , . dr wi 
inh cl, 
lel (7+ 1)2(7+ a) 
ili 0 
i1- ‘ ee ae . 
whence it follows that H, is singular only when ¢% is zero or infinity and is therefore 
ae. reg lear 2/ 
regular in D. Near w& 0 H, = mph 2. 
as 
, 1 
e] | and near 1 /y 0 H, 1 3 8 Y- 
“s | 
Further, on the cut in the s-plane & —x, x > 0, and H, = L+iM, where 
no 
is | 4 x? dr r xi dr 
M ———— and L= >_< 
. | (r+1)'(r— x) (7+1)*(x—7)#" 
ity ere 7X ; x~% 
according as % is on the upper or lower side of the negative real axis of ws. Now if 
as } 8 €+e¢, where &, e are real, |[&| < 1, and € is small, 
lso | ote Qie&c? 
a ~~ 2 2\ 1 72 2)2? 
at a*(1—€?) © a*(1—€?) 
te so that we can correlate positions on the cut in the s-plane with positions on the 
“a negative real axis of & in the &-plane. Thus, if € > 0,¢€ > 0, weare on the upper side 
( . : “ P : 
of the negative real axis of. Let us consider J first; whens is small I ~ —s(2-+-7ci/2a) 
he if res Oand I ~ —s(2—7ci/2a) if res < 0. When s?+1 is small 





(s inf “- (s—i)® log P | 
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and on the cut reJ = +M(1—&) > Oif e€ > 0, and < 0 if e€ < 0. Hence we can 
trace the change in arg J round D. On the infinite circle arg I increases by 27, round 
the small circle at s = —i it decreases by 27, while on the rest of the cut J never 


crosses the real axis of s, so that the net increase in arg I is zero. Now J is regular 
in D and therefore it never vanishes in D. 
Hence also J has no poles in D and is regular there. Further, near s = 0, 


— 77 8C 


J = ——— if res > 0 
78C ; 
and J = —_—— ifres < 0; 
near s?+ ] 0 J = —s(s— if a S(s— i)? log = | : 


and on the cut im J T - 


(2€+(1—) LP + MXl—SP 


according as e€ is positive or negative. Hence arg J increases by 47 on the infinite 


circle, decreases by nearly 27 round the small circle at s = i, but does not cross the 
positive real axis of s there, and crosses this axis twice near s = 0, once on either 


side of the cut. Hence, the net increase in arg J is zero and J has no zeros in D. 
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SUMMARY 
A rarefaction wave is produced in a two-dimensional channel by the impulsive 

retraction of a piston. In a finite section of the channel there is a slight variation 

in cross-section, and the effect of this variation on the rarefaction wave is considered. 

It is found that, after passing through the transition section, a perturbation is super- 
sit imposed on the main wave which tends to decrease the flow velocity by a term which 
the decreases inversely as the time for a point moving with the wave, but increases linearly 
ther with the time at a fixed point in the channel. In the transition section the distur- 
), 


j bances accumulate and tend to increase the flow velocity if there is a throat, even- 
tually producing a sonic régime. An attempt is made to describe the establishment 
of steady subsonic-supersonic flow in a converging-diverging nozzle. 


Introduction 
THERE is a now familiar one-dimensional theory of the steady flow in a 
duct of varying cross-section which, within its limitations, gives an adequate 


en ele 


picture of the flow in, for example, the supersonic wind-tunnel of inter- 
mittent type. This last-mentioned tunnel is usually equipped with a valve 
at the downstream end, which opens into a vacuum reservoir. When the 
valve is opened, a rarefaction wave proceeds through the working section, 
the nozzle, and eventually into the upstream section. There is thus a certain 
period in the initial stages of the motion during which the flow is unsteady, 
; and experimental evidence suggests that this unsteady flow persists until 
a ‘starting’ shock is formed at the throat of the nozzle. This shock eventually 
proceeds downstream leaving a steady-flow régime behind it. 

It is with the initial period of unsteady flow that the present paper is 
concerned. In order to produce a tractable problem the duct is assumed to 
be a two-dimensional channel, the variations in width being small and 
occurring in a transition section of finite extent. Outside the transition 
section the channel is of uniform width. It is also convenient to replace 
the valve by a piston which is impulsively retracted at the ‘escape’ speed 
of the fluid (see below or (1)). This will simulate the effect of the valve 





and vacuum reservoir in the sense that a centred rarefaction wave travels 
upstream, and there will be no reflection at the piston of the subsequent 
disturbances produced in the flow. 


(Quart. Journ. Mech. and Applied Math., Vol. VII, Pt. 2 (1954)] 
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The basic equations and their solution 

Let the origin of the coordinate system O(x,y) be taken at the initial 
position of the piston with the x-axis along the central line of the channel. 
Let the piston be retracted at time ¢ = 0 with the constant impulsive speed 
Co(1—p*)/u?, where cy is the velocity of sound in the undisturbed gas and 
pu is related to the adiabatic index y by the equation p? = (y—1)/(y+1), 
Until the transition section is reached, a one-dimensional unsteady flow is 
produced by the resulting rarefaction wave in which the flow velocity « 
and the velocity of sound c are given (1, p. 104) by 


u = (1—p?)(xt-1—c,) ) (1) 
c= p2xt1+(1—p2)e,) 
It will be seen that at the piston wu = —c,(1—p?)/p?, and the above choice 
of piston speed is just sufficient to reduce the sonic velocity (and incidentally 
the density) to zero there. It is called the escape speed (1) and completes 
the rarefaction wave since the latter then ends in a vacuum. 

In the disturbed flow let V be the velocity, ® the velocity potential, and 
p the density. Then the equations of conservation of mass and energy are 


P 4W.(pV) = °| 
ot 

i ; (2) 
pe 


0 yo ea ao} 
ac °° 2 


and elimination of p from equations (2) yields the relation 


0. (- +V.V)IPe+(V.0) = 720, (3) 
ot? ot ct 


It is not difficult to verify that the velocity potential of the flow produced 
by the rarefaction wave alone is 


$(1—p*){c8 t+ 27t-1— 2c, 2}. (4) 
We thus assume that 
® = }(1—p*){c§t+a7t-1— 2c, 2}+¢, (5) 


where ¢ is the velocity potential of the perturbed flow, and linearize 
equation (3) on the basis of a small perturbation superimposed on the 
main flow produced by the rarefaction wave. After some algebra one 
obtains the following equation for 4, 
py tat-h{ (1 —2pu?)at-!— 2(1—p2*)e9},,—{p2at-1+ (1—p? Jeo}? + 
+2(1—p?)(xt-*§—Co)hy— 2(1—p*)eg tf, + 2p7t-*h, = 0, (6) 
where suffixes denote partial derivatives. 
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Finally, making the substitutions 


é = p?- . (l—p?), ny = wat(l—p*)eot, C = (1—2p?)ty, (7) 
Co 
Fa : 
equation (6) reduces to Pin ~ x en byr (8) 
2u? 
where ¥ = e _ (9) 
1—2p* 


ind suffixes denote partial derivatives. 

Here, € varies between 0 at the piston and 1 at the head of the wave, and 
the corresponding limits for » are 0 and ¢of. 

Since the velocity of sound is zero at the piston, no disturbance can reach 
itand so ¢ and its derivatives must be zero at € = 7 = 0. Also, ¢ and its 
first derivatives must be zero at the head of the wave, where € = 1. Finally, 
¢is an even function of ¢, and if the width of the channel is 2{b+/(z)}, 
where f(a) is uniformly small, then 


Cd 


© = (1—p\(et—ey)f (@) (10) 


when y = 6. In terms of the new variables this becomes 


: 2 2 
r= — tae tor 1-H} (11) 
\p? 
when € = B = (1—2y?)!D. 

The problem is now reduced to mathematical terms. A solution of 
equation (8) is required subject to the appropriate boundary conditions. 
This is most easily obtained with the help of two transformations. The 
variable € is first suppressed from equation (8) by means of a Mellin trans- 
form in which 


1 
b(s) = | (Ege dé, (12) 
0 
, A+ ico 
$6) = 5 | Welerds (13) 
\-ica 


2, p. 46). The path of integration in (13) lies to the right of all the singu- 
lavities of ¥(s). These formulae are easily obtained from the well-known 
relations used in the Laplace transformation by a simple change of variable. 
Since ¢, is zero at both € = 0 and € = 1, equation (8) becomes, in terms of y, 


XS 


Pant ain _ pyr. (14) 
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Lastly, { can be suppressed with the aid of a finite Fourier cosine transform 
(3, p. 78) defined by B 


F(n) = | ¥(L)cos(nmt/p) af, (15) 
$= 5FO+5 > F(n)oos(nzt/8). (16 


Since 7; is zero at ¢ = 0 and prescribed (through equation (11)) when { = 8, 
it follows that equation (14) becomes, in terms of F, 


n? 7 


a. (17) 


OS 1, 
Fn ttn - (—1)"bg— Be 
where yg has been written for (@//0C);_,. 
With nz/8 replaced by k, the general solution of (17), for n + 0, which 
has the correct behaviour at 7 = 0, is 





7 
F= (—1) "gary Ket) | {Sy (as—-1(k) Vy as—(k2) —Yy(as—v(k) Sy (s—1)(k2)} , 


0 x zias+Dyh, dz. (18) 
The interpretation of this expression is difficult. However, in the develop- 
ment of the theory of steady flow in a duct, the assumption is made that 
the flow variables are approximately constant over the cross-section of the 
duct. This is a reasonable assumption if the variations in cross-section 
are small, and the equivalent simplification in the present problem is 
concentration on the first term of equation (16), which gives the average 
disturbance over the channel cross-section. The expression for F(() is 
much simpler, for the solution of (17) with n = 0 is 


v 


U] n 
F(0) = [ y-%8 dy [ u*ib_a(u) du = —— | fg(w)[u%snt-%8 —u] du, 
0 


0 0 (19) 
where account has been taken of the boundary conditions to be satisfied 
at 7 = 0. 

It is convenient first to interpret equation (19) when the channel wall is 


defined by the equation y = b+H(x—a), (20) 


where H(x) is the Heaviside unit function. The results for an arbitrary 
variation in cross-section may then be deduced by superposition. 
The substitution f(~) = H(x—a) in equation (11) gives, for dg, 


— _ to(l—p*) = 5[2(1—--*) —al 21 


y(1— 2p?) 
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nsform | where 5 is the Dirac delta function. Using equation (12), 4g may now he 
determined. Explicitly, it is given by 











(15 c.(1—p2) ; 1— 2 
; alll 5} 7 —e $— 
$e = —— (1—é)8 (1 )—a &s-1 dé 
pe(1—2y*)* J \p g J 
(16 : 
| c,(1 2)s+1,,2/ 9a —8—2 a 
~~ ez, (i u? 1— “\H(y—a). (22) 
. ; n(1—2p*)® 7 n 
c=, | 
| When this is substituted in equation (19), the resulting equation for F'(0) is 
- , 2\s+1,,2 a —s—2 
(17 F(0) J Co(1—p*) ey i fy — MN fgnegt—ao_ gp 
(1—as)(1—2p?)* J | u) | 4a u 
a 
(23) 
which Now, by (13), the interpretation of 7°/(s—1/a) is 
} r , 
| - ((&\-* ds "| ieee 
| x | (J = (2) He—8. (24) 
2m J tT} 8s—l/a T 
A 
(18) It follows that the interpretation of (23), with the help of equation (9), is 
elop- Cy( 1 — po?) +8 9(1 — 2?) [ l—a/u 
that Della (1—pa/u)? ste’ 
rf the a 
tic | — ge? x 1—p* 
cti m H! | u _¢l _ HI! an 108 an) el du. (25) 
m is \1—p2a/u\n J \i1—p2a/u*J 
i i The term independent of y in the Fourier series for ¢ is just (25) multiplied 
Y) 1s by 8-1 or 1/(1 27)*d. 
When 7» x, (25) is identically zero. Interpreted physically this means 
| that there is a region of length a/u in the neighbourhood of the piston which 
is not affected by the disturbances. It is, however, the region x > 0 which 
(19) is of practical interest, and we consider first the case 
1—p? 
sfied —— et -. (26) 
1—p*a 1) 
lis | In terms of our original variables, this condition, by equations (7), is equiva- 
lent to 2 
(20) l we 1—p* ; 
Cot p?(x—a)+(1 p)egt 
ary 
OI x ra. 
Now, fora < u n, the expression 


21) | 1—p* _(u\* 
1—pa u\n 
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is monotonic increasing and lies between 





a —py2 
(* and z. P —. 
v7) 1—p*a/y 
Thus, when (23) is satisfied, 
Al 1—p? u *_ | 
\1—p?a/u\n 
is zero throughout the range of integration, and the main term in ¢ is, from 
(25) multiplied by 1/(1—2y?)%d, 


pa€/{E—(1—p2)} 








d Co(1—p?)i +o 1—a/u du 
7 bg" J {1—p?a/u}?+vo 
1 
= , _ H+(1—p*)cot/x 

C) a(1—p?) i+ . z, - —1/a j—-1 du 

= — ——— | pu“ — -++- | —p* — - : 
2b Cot ; I (1 peu-l)? +1) 
1 
(27) 


For an arbitrary variation f(x) in the width of the channel, we may write 
oO 
je) = f'(a)H(«—a) da, 
0 


and, by superposition, equation (27) can then be generalized to 


2)1+1/« . —1/a 
od = Se nies { af'(a) da(u = + 1—*) x 
. BP+(1—p*)eot/x 

E 1—u-! 
x rennin Ha. TOA 
(1—piu-t)Prue 

At all points beyond the transition section, where the channel is again of 
uniform width, this is the only contribution to the potential. Also, in this 
case, the dependence on z of the first integral occurring in (28) is spurious, 
and the potential is actually a function of x/(c,t). Consequently, ¢, is of the 
form af'(a)da_ (x ‘a 
= ele yi 

(say), and is positive for a nozzle outside which there is no net change in 
channel cross-section. Since the main flow is in the negative x-direction, the 
tendency is for the flow velocity to be reduced by the amount given by 


(29). Furthermore, by the second of equations (2), we have 





c= ¢ oe {(l—p?)(xz/t—c))+¢,,}" ap" {$(1—p?)(c§—x?/t?) +4} 
ee a BN oT Pss —;5 r Bo tS 
(30) 
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To the first order of small quantities, this gives 

p*{x/t—co}b,+{n?/(1—p*)}dy (31) 
pra /t+-(1—p*)eg . 


ind, when ¢ is a function of /(cyt), this simplifies still further to 


Cc fa t+(1 jt*)Co} — 





c = {u2x/t-+(1—p2)c,}+—" — ¢,. (32) 
1—p* 


This gives the increase in sonic velocity corresponding to the decrease in 


flow velocity. 
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For y = 1:4, w? = 4, a = }, and the last integral in (28) can be evaluated 


in terms of elementary functions. The result is not exhibited explicitly, 
for in practice the function G of (29) is the more interesting quantity, and 
its variation is more easily appreciated by means of a graph. Fig. 1 shows 
the variation of G with 2/(cgt). 

The head of the wave corresponds to 2/(cyt) = 1, where G = 0. As the 
wave is penetrated G increases monotonically and satisfies the asymptotic 
equi ? 

juality “in 25 ce 


— > 
i2 z* 


aS x/(Cyt) > 0. At any fixed time, the variation of G in the wave represents 
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the variation in perturbation velocity, apart from a constant factor (see 
equation (29)). When the factor is positive, as in the case of a nozzle, this 
will correspond to a decrease in the magnitude of the main flow velocity, 
which is in the negative x-direction. For a point moving with the main 
rarefaction wave (a/(c)t) = constant), equation (29) shows that the per- 
turbation decays through a factor ¢-!, but the asymptotic behaviour of ¢ 
implies that the magnitude of the perturbation increases linearly with the 
time at a fixed point in the channel. This is because, at any fixed point, 
the flow speed in the main wave is steadily increasing in magnitude to the 
value (1—,”)cy, and the sonic speed is steadily decreasing to the same value 
(see equations (1)). This has the effect of increasing the intensity of the 
perturbations produced by the variations in cross-section (since the flow 
speed increases) but decreasing the speed with which they are propagated 
upstream (since this is the difference between the flow speed and the sonic 
speed). Consequently, there is a tendency for the perturbations to accumu- 
late. This will cease, however, as soon as sonic conditions are attained at 
the throat, for then no disturbances will be propagated upstream, and the 

flow will then become steady, at least in the upstream section. 
We now return to the basic equation (25) and consider the case 
1—p? » 
§&< ——-, (33) 

1—p*a/n 
which is equivalent to x <a. The function H{(1—,?)/(1—p?a/u)—€} is 
now equal to unity throughout the range of integration, and the contribu- 
tion to ¢ (which is (25) multiplied by 1/(1—2y?)!) will now be 


2 


7 : 
C(I =n) af ¢ i | ie (1 —* w) 
2bé1/% ‘ (1—p2a/u)?+Vo 


Uy 
Uy 
+ 


__ 4,2) +o ” 
- __ €o(1—p*) me : u) - (34) 
(1—p2a/u)?+t 





where wu, satisfies the relation 
I—p? /u,\* on 
a (") ax £, (35) 
1—p*a/u,\ 7 
It will suffice to consider the behaviour of (34) for large values of t, which 
implies large values of y. In this case, from (35), 


2 1 
Uuy~ (1 —n*)egt-+0(—) 
0 


for a fixed value of x. Thus, (34) is asymptotically equal to 


— ct) [1 —p)e9t—a(d—p?)log co]. (36) 
<0 
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For an arbitrary variation in the width of the channel, this will generalize to 


« 


_ Col ae (1—p*)egt | f'(a) da —($—p*)log cot [ as"(a) da - (37) 


zx 
The complete asymptotic expression for ¢ is obtained by adding to (37) the 
leading terms in the asymptotic expansion of (28). When this is done, one 
finds that 


r 


b~ - a Ke) (1- 2) oo! af'(a) da —($—p?*)log cot af'(a) da +- 
20 x J 
0 


© 


0 


+-(1—p?)egt f' (a) da —(4—p* log cyt | af"(a) da . (38) 
Differentiation of this expression gives 
d,~ == H')"cot [ af'(a) da. (39) 
2ba® = 


0 
In the region downstream of the variations in cross-section, where 
f(z) = 0, (39) is zero and there are no large-scale disturbances. But in the 
case of a converging-diverging nozzle, with no net change in channel width 
zx 
on the two sides of it, the factor 2-? { af’(a) da will, in the nozzle near 


0 
the downstream end, be negative and decreasing as x increases (for example, 


if f’(2) k for x > c, then the last-mentioned factor is —k(1—c?/x?) 
and is monotonic decreasing). This will tend to increase the flow speed, 
which is in the negative x-direction. However, this tendency must be 
reversed before f’(x) = 0, that is, at some point slightly downstream of the 
throat. As x increases further, ¢, will begin to increase and, since f’(x) 
becomes positive beyond the throat, ¢, will eventually become positive 
and attain a maximum at some point before f’(x) again becomes zero, 
that is before the end of the nozzle is reached. Thereafter, for a given 
value of ¢, ¢, will decrease as x increases. Now when ¢, is positive and 
increasing in the direction of positive x, the perturbation on the flow speed 
in the main wave is negative and increasing in the direction of flow. Thus, 
as the nozzle is entered from the upstream side, the perturbation on the 
flow speed is negative. It decreases, reaches a minimum, increases, and 
becomes positive before the throat is reached. Maximum speed is reached 
just downstream of the throat, and beyond this point the perturbation 
remains positive but becomes smaller (for fixed ¢) as the downstream end 


of the nozzle is approached. 
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The inference is that a sonic speed first appears somewhat downstream 
of the throat. There will most likely be further adjustment of the flow in 
this region due to non-linear terms which will become prominent as the 
disturbances accumulate. Once a sonic régime is established, however, 
there will be no further disturbances propagated into the upstream section 
of the channel, and the linear increase with time of the perturbation 
velocities in this section will cease. The flow in the upstream section will 
then consist of a rarefaction wave plus perturbation effects (for a point 
moving with the rarefaction wave 2/t is constant and the perturbations 
will decay through a factor ¢-1, see equation (29)). At the tail of the 
rarefaction wave will be a region of steady flow which grows in time and 
extends upstream from the nozzle. But the flow in the downstream section 
is by no means steady, and further effects occur which are outside the scope 
of the present theory. When sonic velocity is established near the throat, 
the downstream section of the nozzle is acting like a diffuser, and diffuser 
flow is unstable. Experimental evidence (4) shows that a shock is formed 
just downstream of the throat. The combined effect of the shock and the 
adverse pressure gradient in the diverging section of the nozzle causes the 
flow to separate, the separation beginning at the junction of the wall and 
the shock. As the pressure falls, the shock moves downstream taking with 
it the points of separation (at each wall of the channel) and leaves behind 
it a steady supersonic régime in the downstream section. 
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